Lattice Gas Automata for Reactive Systems 



Jean Pierre Boon 
Center for Nonlinear Phenomena and Complex Systems 
Universite Libre de Bruxelles - Campus Plaine - CP. 231 
1050 - Bruxelles, Belgium 

David Dab 
Department of Chemistry 
Massachusetts Institute of Technology 
Cambridge, Mass. 02139, USA 

Raymond Kapral 
Chemical Physics Theory Group 
Department of Chemistry 
University of Toronto 
Toronto, Ontario M5S lAl, Canada 

Anna Lawniczak 
Department of Mathematics and Statistics 
Guelph- Waterloo Program for Graduate Work in Physics 
University of Guelph 
Guelph, Ontario NIG 2W1, Canada 
(October 1995) 



(to appear m PHYSICS REPORTS) 



Reactive lattice gas automata provide a microscopic approach to the dynamics of 
spatially-distributed reacting systems. An important virtue of this approach is that it 
offers a method for the investigation of reactive systems at a mesoscopic level that goes be- 
yond phenomenological reaction-diffusion equations. After introducing the subject within 
the wider framework of lattice gas automata (LGA) as a microscopic approach to the phe- 
nomenology of macroscopic systems, we describe the reactive LGA in terms of a simple 
physical picture to show how an automaton can be constructed to capture the essentials 
of a reactive molecular dynamics scheme. The statistical mechanical theory of the au- 
tomaton is then developed for diffusive transport and for reactive processes, and a general 
algorithm is presented for reactive LGA. The method is illustrated by considering appli- 
cations to bistable and excitable media, oscillatory behavior in reactive systems, chemical 
chaos and pattern formation triggered by Turing bifurcations. The reactive lattice gas 
scheme is contrasted with related cellular automaton methods and the paper concludes 
with a discussion of future perspectives. 
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I. INTRODUCTION 



Often systems with a large number of degrees of freedom exhibit macroscopic behavior where the 
details of the microscopic dynamics are relatively unimportant. This feature makes a general description 
possible: systems with different microscopic characteristics can be described on the macroscopic scale 
with a generic set of equations where the specific nature of the components of the system is refiected in 
the numerical values of a restricted number of coefficients which enter through the constitutive relations. 
When the macroscopic equations are recast into a non-dimensional form the global behavior of the 
system depends on a limited number of universal control parameters where the microscopic nature of the 
elementary constituents does not appear explicitly. The macroscopic level of analysis therefore provides 
a phenomenological description where (i) the complexity of describing the dynamics of the microscopic 
degrees of freedom is bypassed, and (ii) many different physical phenomena are grouped into a limited 
number of classes. Classical fiuid mechanics offers a striking example through Reynolds' dynamical 
similarity law (1883). 

Nevertheless, it is obvious that in this process of reasoning the connection between the phenomenol- 
ogy and the underlying microscopic mechanisms has been lost. Consequently, the macroscopic level of 
description, which employs average quantities, cannot be used to analyze how large scale phenomena are 
triggered by local and/or transient deviations from these averages: it neglects fiuctuations. This is one 
of the objects of the statistical mechanical approach which establishes the microscopic basis of the prop- 
erties of many-body systems. We have observed that it is not necessary to have a complete knowledge 
of the details of the microscopic interactions to understand how macroscopic phenomena emerge. We 
shall demonstrate that it is possible to develop a mesoscopic approach to macroscopic phenomena by 
modeling the microscopic dynamics by means of a simplified description, provided the basic requirements 
of fundamental physics and chemistry are correctly incorporated, i.e. the conservation laws, symmetry 
properties, reaction mechanisms, etc. 

This philosophy has been exploited successfully in hydrodynamics and in statistical hydrodynamics 
where schematic microscopic models - lattice gas automata (LGA) - were initially developed with the 
goal to produce simplified models for statistical mechanics [1] and subsequently to provide new com- 
putational tools for the study of complicated problems in fiuid dynamics [2]. Various LGA have been 
proposed to model macroscopic physical phenomena such as 2-d and 3-d fiows at moderate Reynolds 
number, immiscible multi-component fiuids and fiows in porous media [3]. These LGA share a common 
structure where point-like particles move along the links of a regular lattice and interact on the nodes 
through collisions conserving particle number (mass), energy and momentum. The key point is that 
these conservation laws along with some symmetry requirements for the lattice structure suffice as basic 
ingredients for the emergence of macroscopic fiuid behavior in agreement with Navier-Stokes equations.^ 

Quite naturally attempts were also made to construct a "lattice chemistry" in very much the same 
way that hydrodynamics had been successfully implemented on lattices. The study of the combination of 
hydrodynamics and chemistry via LGA methods could eventually lead to a simplified approach to highly 
complex problems such as reactive fiows and combustion [4]. The operational feasibility of simulating 
reactive fiows was demonstrated by implementing a simple specific reactive scheme (A + B^2C)ma 
two-dimensional lattice gas subject to a shear constraint [5]. 

The goal of the research described here is different in scope. We present a general LGA approach 
proceeding from micro-dynamical equations governing the dynamics of reactive "particles" to the macro- 
scopic behavior of reaction-diffusion (RD) systems. An important distinction between macroscopic fiuid 
dynamics as described by the Navier-Stokes equations and RD systems is that for the latter there is no 
unique similarity law. As a consequence, universality cannot be achieved in the same way for reactive 
phenomena. The LGA scheme for reactive systems therefore has as its goal the construction of a mi- 
crodynamics for a class of reaction-diffusion systems which exhibit space- and time-dependent solutions 



^Galilean invariance is usually not satisfied at the microscopic level of LGA and this affects the macroscopic 
behavior. However, Galilean invariant macroscopic behavior can usually be obtained through scaling for quasi- 
incompressible, single-species flow [2]. 
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corresponding to a broad variety of observed and predicted phenomena in reactive systems. Our appli- 
cations include bistability, chemical wave and pattern formation processes (for example spiral waves and 
Turing patterns) and systems that exhibit chemical chaos. 

Our approach utilizes the concepts of statistical physics which provide a mean to establish the connec- 
tion between microscopic and macroscopic properties of many-body systems. In this respect, an important 
feature of the LGA approach is the ability to yield a mesoscopic level of description of the space- and 
time-dynamics which can be used to investigate the role of fluctuations in spatially-distributed reacting 
systems. The analysis of fluctuation correlations appears to be crucial in RD systems when macroscopic 
phenomena are triggered by the amplification of microscopic fiuctuations. Such fiuctuations are intrinsic 
to LGA because particles have a discrete nature in this approach and because the rules governing the 
dynamics are probabilistic. For instance, it is possible that global solutions obtained from the automaton 
may be at variance with the macroscopic behavior predicted by the phenomenological equations. The 
origin of such discrepancies stems from the fact that phenomenological descriptions do not incorporate 
naturally molecular chaos effects. One is then led to analyze the limit of validity of classical macroscopic 
reactive kinetics and consequently to address fundamental questions underlying the statistical mechanics 
of reaction-diffusion systems. The simplicity of the microscopic dynamics of the automaton also yields 
operational advantages with respect to classical approaches using fioating point algorithms. From an 
operational view point, LGA provide "stable" algorithmic prescriptions for the simulation of the class 
of reactive phenomena considered^ . Furthermore the LGA method offers interesting perspectives for the 
investigation of complex macroscopic phenomena which are difficult to treat with classical analytical or 
computational methods and for situations (a.o. complex media, complicated boundaries, etc.) where 
quantitative laboratory experiments are difficult. 

Of course, other methods have been used to explore the dynamics of spatially-distributed reactive 
systems from a more microscopic perspective than that provided by macroscopic field equations. Full 
molecular dynamics gives the most detailed classical description of the reactive dynamics. Simulation 
of large systems, where macroscopic spatial or temporal structure is possible under far-from-equilibrium 
conditions, using realistic molecular potentials for the scattering events, is still beyond the scope of existing 
computational power. However, for model hard sphere systems where the reactive events are basically 
"coloring" processes there has been some work on the effects of fiuctuations on limit cycle oscillations 
and other aspects of bifurcations in far-from-equilibrium conditions. [6-8] While these simulations are 
limited to relatively small numbers of particles compared to the automaton simulations, they have the 
advantage that the energetics of the reactive events are treated in a more sophisticated manner. 

The modeling strategy and perspective taken here are also similar to those for random walk model of 
reacting systems. These random walk models have been used in a number of studies to explore specific 
reactive systems and the validity of mean field descriptions. [9,10] This class of models is contained within 
the general framework presented here. 

In addition, reactive lattice-gas automaton models have a number of features in common with birth- 
death master equation models for reactive systems. [11-14] Birth-death master equation descriptions 
have been and continue to be used to gain an understanding of the role of fiuctuations on far-from- 
equilibrium rate processes and should be viewed as an approach which is complementary to the lattice-gas 
automaton methods described here. One may also compare reactive lattice-gas automata to kinetic Ising 
models where one has discrete space-time dynamics with discrete spin variables. [15] Such models, in very 
simplified contexts, can be used to discribe reaction-diffusion dynamics at a mesoscopic level. 

The paper is organized as follows. We first present in Sec. II the physical picture behind the automaton 
rule construction and an overview of the types of results one may obtain through its use. 

The next two sections give a statistical mechanical description of the automaton for a simple one- 
variable system whose dynamics occurs on a square lattice. This allows us to introduce and illustrate 
a number of features involved in the model construction and use. Section III on diffusion is devoted 
to the description of non-reactive dynamics where the passage to a discrete Boltzmann equation, and 
subsequently to the diffusion equation, is made. Section IV considers the modifications that occur when 



^By stable we mean that the algorithms are exempt of numerical overflow and roundoff errors. 
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reactive processes also take place. The combination of diffusion and reaction yields the full automaton 
dynamics. In this section we describe the passage via the lattice Boltzmann equation to the linearized 
reaction-diffusion equations and the compatibility with the phenomenological description. 

The general reactive lattice-gas automaton rule is stated in Sec. V. After discussing how the consid- 
erations of the previous sections can be extended to the multi-species case, a detailed strategy for the 
construction of the reaction probability matrix, which forms a central part of the rule construction, is 
given. The combination of this implementation of the reactive step with variants of the diffusion dy- 
namics allows a variety of systems to be treated. This section constitutes a "receipe" for algorithmic 
implementations of the automaton. 

Specific applications are treated in sections VI to IX where we present some of the most typical examples 
of reaction-diffusion phenomena in bistable and excitable media, Turing bifurcations and chemical chaos. 
The next section. Sec. X, briefiy considers the generalization of the model to cases without exclusion, 
which allows a broader class of systems to be described. Section XI briefiy comments on the relation 
of the reactive lattice-gas automaton considered here to other cellular automaton methods for reactive 
systems. The final section of the paper gives some perspectives for future work in the field. 

Hitchhikers guide to the contents: This article was written with several types of reader in mind; so it is 
perhaps useful to provide a road map to its contents. For those who wish merely to gain some familarity 
with the ideas behind the method and its applications to specific problems we suggest reading Sees. II 
and VI to IX. For those who want to implement the method on a computer for their own purposes 
we recommend Sees. II and V-X. Readers interested in the statistical mechanics of the method should 
consult Sees. II-IV. Naturally, dedicated scientists are invited to read the entire paper. 
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II. PHYSICAL PICTURE 



Consider a reactive system where various species X, A, B, C, ... diffuse in a soivent S. Suppose that the 
species X undergoes one or a series of reactions of the form 

aX + ... ^ (3X + ... , (1) 

where the ". . ." represent terms involving the A, B,C, ... species. We assume that these A, B,C, ... 
species have their concentrations kept at constant values through external constraints; thus, the system 
is maintained in a non-equilibrium state. Then the macroscopic dynamics of species X can be described 
by a reaction-diffusion equation 

^^^^ = F{px{v,t)) + DV'px{v,t) , (2) 

where px{v,t) is the density of species X at point r at time t. In this description, the species A,B,C, . . . 
do not appear explicitly: their concentrations are incorporated in the reactive rate F(px). The solvent, 
whose microscopic role is to scatter molecules through elastic collisions, is also hidden in this description 
where it manifests itself only through the diffusion term DV'^px. Now we show that a mesoscopic 
approach which ignores the dynamics of A, B, C, ... and S molecules and retains solely the dynamics of 
X molecules can be constructed to yield the same macroscopic phenomenology as described in (2). 

We consider a rf-dimensional lattice with coordination number where molecules are modeled as point 
particles which undergo displacements along the links connecting nodes. Particles move on the lattice 
with discrete velocities, that is they hop at discrete time steps from a node to one of the neighboring 
nodes as dictated by the particle velocities. Each node of the lattice possesses m channels where particles 
can reside, and an exclusion principle forbids more than one particle to reside in any channel. Therefore, 
the total number of particle at a node is restricted between and m. A velocity is associated to each 
channel so that a particle has the velocity of the channel on which it resides. In most of the models 
we consider in this paper, the number of channels m is equal to the coordination numbers of the lattice 
rric and the allowed velocities correspond to the jumps from one node to a nearest neighbor node in one 
discrete time step. 

The time evolution of this single species mesoscopic system occurs at discrete time steps and follows 
from the iterated application of an evolution operator £ (also called the rule of the automaton) 

[State at time {k -|- 1)] = 5 [State at time k\ (3) 

which can be conveniently decomposed into three basic operations: propagation P , velocity randomization 
or mixing R , and chemical transformation C [16] 

£ = CoRoP . (4) 

These operations will now be described in their simplest forms; generalizations are given in Sec. V. 

(i) During propagation, each particle hops from its channel to the corresponding channel of a neighbor 
node as dictated by the particle velocity; it is a free streaming process in which the number of particles 
and their momenta are conserved. 

(ii) In the velocity randomization step, the velocity configuration is randomly shuffled at each node 
where X particles are redistributed amongst the channels. This operation conserves the number of 
particles at each node but the momentum is not conserved. The momentum changes can be viewed as 
elastic collisions between X particles and a solvent not described on the lattice (a stochastic momentum 
reservoir). In this sense, the role of the velocity randomization is to modify the velocity distribution in 
much the same way as repeated collisions with solvent molecules would do. 

(iii) During the chemical transformation, X particles are created or annihilated at each node in reactions 
of the form 

aX PX . (5) 
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This operation is performed independently and simultaneously at each node of the lattice where a configu- 
ration with a X particles is transformed into a configuration with [3 X particles, with probability P{a, [3). 
Except for a = [3 (which leaves the state of the node unchanged), the reaction operator conserves neither 
particle number nor local momentum. 

An illustration of a sequence Co Ro P of the automaton microscopic dynamics is shown in Fig. 1. 



FIG. 1. Illustration of the microscopic dynamics of a reactive lattice gas automaton on a square lattice. Arrows 
indicate the presence of particles with their corresponding velocity vectors. The figure shows the successive 
transformations corresponding to the evolution operator £ = C o R o P: (a) state at some discrete time k; 
(b) state after propagation P; (c) configuration change after velocity randomization R; (d) state after chemical 
transformation C: this is the state at discrete time k + 1. Boundary conditions are periodic. 

Assume, for a moment, that no reaction occurs (i.e. the reactive collision leaves the system unchanged: 
C = identity). Then, through the repeated application of the operator Ro P^ X particles execute random 
walks; it will be shown in Sec. Ill that diffusive behavior follows in the macroscopic limit where the 
system dynamics is well modeled by the diffusion equation 

^^ = i^vV(M). (6) 

Now, when reactions take place, the number of particles is not conserved and a source term enters the 
diffusion equation to yield (2). If reactions are not frequent (_P(a,/3) ^ 1 for a ^ /3), we shall see in 
Sec. IV that the source term takes the form of a reactive rate 

F{px) = ficx) = Y.{p-a) M (1 - cxT- P{a, /3) , (7) 

where cx = px/rn is the A-particle density per channel (the channel occupation probability). The 
physical interpretation of (7) is as follows. When reactive processes are infrequent, particles on the 
lattice can be assumed to be distributed according to an equilibrium distribution. Under this local 
equilibrium assumption, channels are independentely populated, and the probability to have a particles 
(and m—a empty channels) on a node is given by a binomial distribution (™) (1 — cx)™"". Multiplying 
this expression by P{a,j3) gives the probability for producing [3 X particles out of a A particles. By 
averaging the variation of the particle number per node (/3 — a) one obtains the polynomial rate given 

by (7). 
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Generally, different probabilities P{a,j3) produce different reactive rates F{px), which allows one to 
model various macroscopic phenomena. Whether a particular reactive rate F(px) can be realized by an 
appropriate choice of values for the probabilities P(a,l3), and how this choice is set, will be discussed in 
Sec. V. There are restrictions imposed on the types of systems that can be treated with LGA because of 
the exclusion principle, but as discussed in Sees. X and XI, LGA and CA models without exclusion can 
be construted to remove some of these restrictions. 

We close this section with two illustrations of LGA model results that presage the more extensive 
developments in Sees. VL IX. The first example concerns a RD system based on the reaction scheme 




2X + 5 4 3X , (8) 

known as the Schlogl model [17], one of the best-known examples of a reactive system that gives rise to 
bistable states. The model is usualy considered under conditions where the concentrations of A and B 
are maintained at constant values by an appropriate external feed of chemicals. In this situation, the 
reactive rate for X is a cubic polynomial 

F(px) = ao + aipx + a-jPx + ^^aPx y (9) 

whose coefficients Gi depend parametricaly on the concentrations of A and B and on the kinetic constants 
k±i, k±2. For appropriate values of the parameters, the rate equation px = F(px) shows two stable and 
one unstable steady solutions. Figure 2 concerns this bistable regime and shows how a system initially 
prepared in the homogeneous unstable state evolves through domain formation and front propagation. 



FIG. 2. Evolution from the unstable state for the Schlogl lattice-gas automaton. 

Richer varieties of space- and time-dependent phenomena are seen in multi-species systems. Such 
systems are described by sets of RD equations for which automaton models can be constructed. As an 
example we show in Fig. 3 spiral wave pair formation as obtained in a mesoscopic simulation of the 
two-species Selkov model. Lattice gas automaton models for such multi-species systems provide tools for 
the study of the effects of fiuctuations on a variety of chemical pattern formation processes. [18,19] 
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FIG. 3. Spiral wave evolution in the Selkov lattice-gas automaton. 



In the next two sections we pursue the development of the reactive lattice-gas automaton for a single 
reactive species in order to illustrate the passage to discrete Boltzmann and linearized (reaction)-diffusion 
equations. The generalization to multi-species systems is outlined in Sec. V. 

Note: For the sake of keeping with conventional notations, we use the symbol P for the propagator 
operator and also for the probability matrices, e.g. P{a,j3). However no confusion should result since all 
probability matrices are distinguished by their arguments. 
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III. DIFFUSION 



In a non-reactive cellular automaton (C = identity), each particle executes a discrete random walk 
through the repeated application of the propagation and velocity randomization operators. Since the 
random walks performed by different particles are only weakly dependent^ we can expect the macroscopic 
behavior of the particle density to be diffusive. Indeed, this can be formalized and demonstrated as will 
be shown in this section. The approach starts from the microdynamics which is expressed in terms of 
Boolean random variables. An exact equation is derived for the average particle populations and we show 
how this leads to a diffusion equation in the macroscopic limit. Fluctuations and small scale dynamics 
are also considered. For the sake of simplicity, the discussion is restricted to a single chemical species 
whose dynamics takes place on a square lattice with use of a particular velocity randomization operator. 



A. Space, time and system states 

It is convenient to measure the (discrete) time k by the number of applications of the evolution operator 
£. With this time unit, the time interval between two successive states of the system is equal to one and 
k takes integer values; without loss of generality we shall assume that k = labels the initial state of the 
system. 

The physical space C considered here is a rectangular subset of a 2-d square lattice. The rectangle has 
Li nodes along one of its edges and L2 along the other edge (edges are assumed to be parallel to the 
main directions of the lattice). For finite values of _Li and/or L2 we assume periodic boundary conditions 
along the corresponding directions. To label the nodes of £, we introduce a vector r which is conveniently 
projected on a basis set made of two orthogonal vectors connecting a node to two of its first neighbors. 
In such a basis, the components of r = (ri, r2) have integer values. When the system is finite these values 
are considered modulo Li and L2 respectively. Each node of the system can be occupied by particles 
differing only by their velocities which are restricted to the four unit vectors connecting a node to its 
nearest neighbors: {cj- : i = 1, . . .4}. 

An exclusion principle forbids more than one particle with a given velocity to be at the same node at 
the same time. The state of a node is therefore completely defined by four Boolean occupation variables 
whose values specify the states of the four channels, w(r, i), i = 1, . . . , 4, at the node: 

1 if at time k there is a particle 

with velocity Cj- at node r and 
time k; 

rji(r,k)={ (10) 
) if at time k there is no particle 

with velocity Cj- at node r and 
time k. 

When there is a particle with velocity Cj- at node r, the channel w(r, i) is said to be occupied (rii(r) = 1) 
and unoccupied otherwise (rii(r) = 0). 

The four occupation variables {rii(r, k) : i = 1, . . . , 4} which determine the configuration of a node 
will often be denoted as a Boolean vector (or word) 

v{TC,k) = (ryi(r,A;))i=i_^^^_4. (11) 

Here and in the sequel the angle brackets (• • are used to denote the elements of a vector. This vector 
takes its values in the state space S of all 2'*(= 16) four bit words. The state of the full system at a given 
time k is determined by the knowledge of all the occupation variables at that time: 



•^When two or more particles reside on a node, their velocities cannot be randomized independently because of 
the exclusion principle: if the velocity changes were independent, different particles would sometimes end up with 
the same velocity. 
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'n{-,k) = {'n{r,k) ; r e £} . 
The phase space in which this "Boolean field" takes values is denoted by F. 



(12) 



B. Propagation 

In the propagation step, each particle moves in the direction of its velocity vector to the first neighboring 
node where it occupies the channel with the same velocity label. The resulting state transformation 
is equivalent to that produced by a free fiight of each particle during a unit time step. However we 
emphasize that we do not associate any duration to the propagation operation which is simply defined as 
a state transformation; this allows us to use more than one propagation operator to define the evolution 
operator £. 

In mathematical terms the propagation operation is defined by its action on the state variable s(-): 

P : s(.) ^ s^(.) : s,(r) ^ sf (r) = s,(r - c,); (13) 
it is a permutation of the Boolean occupation variables on the lattice C and its inverse is given by 
p-i : s(.)^s^"\.) : s,(r)^sr\r) = s,(r + c,) . (14) 
Similar definitions can be written for the microdynamical variables rji. 

C. Velocity randomization 

A general velocity randomization operator performs a random permutation of the channel occupation 
variables at each node. Here we consider particular velocity randomization operators where the permu- 
tations are selected independentely at each node by a stochastic rule which ignores the node state and 
its location on the lattice. In addition, we require the selection rule to be invariant under the "node 
symmetry group" (i.e. the subgroup of channel permutations corresponding to rotations and refiections 
of the velocities in physical space). The 2-d square lattice model considered here has four channels at 
each node and the number of different ways to permute them is 4!(= 24). Let H be this set of permu- 
tations. The velocity randomization operator is fully characterized when a probability is assigned to 
each possible channel permutation in H with the requirement that Pw = 1, and = when tti 
and 172 are equivalent to within a node symmetry. To express the velocity randomization operator R in 
mathematical terms, it is useful to introduce a set of random variables. Let ^ = (^7r)7rGn be a random 
Boolean vector such that: 

1. for each realization of ^ there is one and only one component tt (the randomly selected permutation) 
such that = 1 ; 

2. for each possible permutation tt the probability of the event = 1 is given by . 

Each time we need to perform a velocity randomization at a node we can draw a Boolean vector ^ 
distributed as ^ and select the permutation tt such that = 1; all selections of ^ are assumed to be 
mutually independent. Suppose that there is only one instance of R entering into the evolution operator 
£ as in (4). Then we need one copy of ^ at each node r and at each time step k. Let ^(r, k) be that copy. 
With these notations, the transformation of the state resulting from velocity randomization at time k 
can be written as 

■sf(r)=Y,Ur,k)Y,Pj^{^>j(r), (15) 
Tren j 

where Si and sf^ denote the pre- and post- velocity randomization states, respectively, and where {pijiir) : 
i,j = 1, . . .4} is a Boolean matrix coding the channel permutation tt in the following way: 
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{1 if the permutation tt maps channel i on channel j ; 

(16) 
otherwise. 

A simple example of such a velocity randomization operator is constructed as follows: the four possible 
rotations (0, 7r/2, tt, 37r/2) are equally probable'* and all the particles residing on a given node undergo 

with probability 1/4 
Hr,k)= { J.'^'-JJ i,{ ^f^^^ (17) 

tbtd. 

independently at each node r and at each time step. We term this the J^-eqm-rotatton velocity random- 
ization. An expression similar to (15) can be written for microdynamical variables. 

Velocity randomization can also be defined in terms of transition probability matrices. Given a velocity 
randomization operator R, we can evaluate for each pair of configurations s and s' the probability 
R(s,s') that s is mapped onto s' by R; these probabilities are collected in a matrix R with elements 
{R(s,s') : s,s' G S}, the (local) transition probability matrix of R. Note that there are transition 
probability matrices which do not correspond to any velocity randomization: to represent an acceptable 
velocity randomization, the matrix R(s, s') must be compatible with (i) the particle number conservation, 
(ii) the "state independent" nature of the permutation selection, and (iii) the node symmetry invariance. 
When the transition probability matrix R(s, s') is obtained we define a matrix whose elements are Boolean 
random variables: 




fes' : eS} (18) 

such that: 

1. for each configuration s there is one and only one configuration s' (the selected post-randomization 
configuration) such that ^ss' = 1, 

2. the probability of the event ^ss' = 1 is given by R(s, s') . 

Each time we need to perform a velocity randomization, we can draw a random matrix [^ss'] and replace 
the node configuration s by the configuration s' such that ^ss' = 1 (all selections are assumed to be 
mutually independent). 

If only one velocity randomization is performed at each time step, then a single copy of [^ss'] is needed at 
each node r and time k; we denote that copy by [^ss'(^) ^)]- With this notation the velocity randomization 
at time k transforms a system state s into s"^ given by 

sfir) = J2U',ir,k)s'/ , (19) 

s" 

which we can also write as 

4 

.f (r) = J2 €s's"(r, k) s'l n(s, (r, k)f^^'^^{l - (r, k)f^-^'^^ , (20) 
s',s" i=i 

where we make use of the convention that = 1 when s = s' {s and s' are Boolean variables)®. 



*The irrational number tt here should not be confused with the same symbol used earlier for a permutation. 
^We recall that the channel indices are defined modulo 4. 
^In this expression the product 
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D. Dynamics and microdynamical equations 



The simplest evolution rule we can construct with the propagation P and the velocity randomization 
R operators is one in which we apply R and P sequentially to make the system evolve from k to k + I. 
This can be done in two different ways: 

£ = RoP, (22) 

or 

£ = PoR. (23) 

The corresponding dynamics 

7i{-,k+X) = RoP-qi-^k) , (24) 

and 

T,(.,fc+ 1) = Poi?T7(-,fc) , (25) 

are referred to as RP— and _Pi?-dynamics, respectively. They are similar but not strictly equivalent 
because R and P do not commute. The main difference is that the system is always observed after the 
propagation step in the _Pi?-dynamics, and after the velocity randomization in the i?_P— dynamics ^. 

Explicit microdynamical equations are obtained by substitution of the analogs of (13) and (15) for 
microdynamical variables for P and R in (24) and (25) yielding 

rii{v, k + l) = ^ ^^(r, k) ^Pii(7r) rij{v - cj , k) , (26) 

ttGII j 

for the i?_P— dynamics, and 

T]i{r + ci,k+l)= ^ ^^(r, k) ^Pji(7r) T]j{r, k) , (27) 
Tren i 

for the _Pi?— dynamics. Alternative expressions follow from using (20) for the randomization, i.e. 

4 

ry,(r, k + I) = ^6s'(r, k) s[ \{ (r - c, , k) (1 - ry,(r - c,,k)f-^^^ , (28) 

and 

4 

ry,(r + c„ k + 1) = ^6s'(r, k) s[ ^'i^, k){l - ry,(r, fc))(i-^^) . (29) 

s,s' j = l 

for RP— and _Pi?— dynamics respectively. 

Since both dynamics are similar, we will restrict the detailed analysis to the _Pi?— dynamics; the 
development for i?_P— dynamics follows exactly the same lines. 



X\{s,t'''(l-s,f-'^^ (21) 
j = i 

indicates whether the equality s = s' is true (value 1) or not (value 0) and therefore selects only the term s = s' 
in the sum over s'. Other expressions can be used as equivalent indicators of (s = s') but it will be seen that this 
particular form is the appropriate choice to derive a Boltzmann equation for reactive models. 
'^It is clear that both dynamics are very similar when the sequence [RoP)" is cast into the form Ro(PoR)^"~^^ oP . 
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E. Probabilistic approach 



Since the random vectors ^ = (^7r)7rGn appearing in the microdynamical equations (26)-(27) are 
(i) equally distributed, (ii) independent of the past evolution of the system and (iii) mutually inde- 
pendent, it follows that the entire evolution process {ri(-,k) : A; = 0, 1, 2 . . .} is a stationary Markov 
chain defined in the phase space F. This process is fully characterized by its transition probability matrix 
with elements 

5(s(.), s'(.)) = r {v{-, k + l) = s'(.)|r,(-, k) = s(.)) , (30) 
= P(5s(.) = s'(.)) , (31) 

which gives for each pair of configurations s(-) and s'(-) the probability to find the system in a state s'(-) 
at time k + 1 when the state at the previous time k was s(-). The Chapman-Kolmogorov equation of the 
Markov chain 

v(s'(-),k+i)= J2 m-),^'(-))m-),k), (32) 

with 

P(s(.),fc) = P(r,(.,fc) = s(.)) , (33) 

governs the evolution of probability measures defined in the phase space and can be viewed as a Liouville 
equation for the lattice gas automaton. Note that in statistical mechanics the Liouville equation describes 
the deterministic evolution of a statistical ensemble of initial configurations; in the LGA models considered 
here, randomness is also present in the dynamics. 

The transition probability matrix £(s(-), s'(-)) of the full dynamical process can be expressed in terms 
of the transition probability matrix of the operators P and R. Since the propagation _P is a deterministic 
permutation of the channel occupation variables of the entire lattice, its transition probability matrix 

P(s(.),s'(.)) = P(i^s(.) = s'(.)) , (34) 

takes a simple form in which each row and each column has one and only one non-zero (= 1) element: 

ri if s'.ir) = Siir - a) Vre£, 
^(K-), «'(•))= ^ Vi=l,...4; (35) 

I otherwise. 

Note that in the argument of V in (34) P is the propagation operator and should be distinguished from 
the transition probability matrix on the left hand side. 

The randomization operator i? is a product of mutually independent local mixing operators acting on 
single nodes. Each element of its transition probability matrix 

i?(s(.),s'(.))=P(i?s(-) = s'(.)) (36) 
can therefore be expressed as a product of local transition probabilities R(s,s'): 

^^(K-),«'(-))=n^(^W'^'w)- (37) 

The transition probability matrix of a particular evolution operator £ defined as a product of R and 
P operators can be obtained from the matrices i?(s(-), s'(-)) and _P(s(-), s'(-)) by taking their matrix 
product. For instance, we have 

5(s(.), s'(.)) = J2 ^"(•))^(«"(-), «'(•)) (38) 

s"(0 

= n^(Kr),s'^-\r)), (39) 
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for _Pi?— dynamics. If we substitute this expression for £(s(-), s'(-)) in the Chapman-Kolmogorov equation 
of the Markov chain we obtain explicit equations for the evolution of a probability measure defined in F: 



s( )Gr rec 



F. Mean values 

Physically interesting quantities are defined as the expectations of functions of the Boolean variables 
rii{v,k). Among them, the average number of particles per channel {channel density) 

#,(r,fc) = E[ry,(r,fc)] , (41) 

plays an important role. Note that because of the Boolean nature of the occupation variable rii{v, k) we 
can identify Ni(r, k) with the probability to find a particle with velocity Cj- at the node r at time k ; the 
density per channel is therefore equivalent to the reduced one body distribution function in statistical 
mechanics. 

Other important average quantities are: the local mass density ® obtained by summing Ni(r,k) over 
the different channels: 

4 4 

p{r, k) = E[^ry,(r, k)] = ^#,(r, k) , (42) 

8=1 8=1 

and the local mass current density 

4 4 

j(r,fc) = E[^ry,(r,fc)c,] =^#,(r,fc)c, . (43) 

8=1 8=1 

It is also possible to define a local mean velocity by dividing the mass current by the mass density 

note that this quantity is not the expectation of a simple function of the occupation variables. 

More involved average quantities can be defined by considering the expectation of a function involving 
occupation variables at different positions and times. Among them, the space- and time-dependent 
density fluctuation correlation function 

G(r, k- r', fc') = E [ (p(r, fc) - E [p(r, k)]) (p(r', k') - E [p(r, k)]) ] (45) 
= E [ (p(r, k) - p{r, kj) {p{r', k') - p{r, kj) ] , (46) 

where p(r, k) = "Y^l-irjii^ , k), is most important and will be considered in detail in the applications (see 
e.g. Sec. VI). 

It is sometimes useful to consider the densities Ni{v, k) as the components of a vector in the Euclidean 
space Eq 

N(r,fc) = (#i(r,fc))fci,^^^4 . (47) 



*In the systems considered here the particle mass plays no role and the mass may be taken to be unity so that 
the mass density is equal to the number density. 
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In this space, the mass density is interpreted as the scalar product of N(r, k) with the particular vector 
whose column form is 



(48) 



which we call mass vector. Similarly the two components ji(r, k) and ^2(1", k) of the current density are 
the scalar products of N(r, k) with the particular vectors 



N,- 



/ 1 



-1 

V 



N,- 




(49) 



which we call the momentum vectors m the directions 1 and 2 respectively. The set N^, Nj^ and Nj^ 
can be completed with a fourth vector 




(50) 



to obtain a complete orthogonal basis of Eq. The projection of N(r, k) along gives the mean value 
of the difference in particle number between the pairs of channels (w(r, 1), w(r, 3)) and (w(r, 2), w(r, 4)). 
We call this quantity the q-density pq{v,k): 



p,{r, k) = E[^(-1)'+Si(r, k)] = #i(r, k) - N^iv, k) + Ns{v, k) - #4(r, k) 



(51) 



G. Lattice Boltzmann equations 



As is usual in statistical mechanics, one faces in LGA theory the complexity of systems with many 
degrees of freedom. Use of the Chapman-Kolmogorov equation (32) recursively allows one, at least in 
principle, to determine the evolution of any initial probability measure defined on F. In practice this is a 
formidable task for any system with more than a few nodes. However, full knowledge of the information 
contained in the probability measure is generally neither necessary nor interesting. The knowledge of 
the one body reduced distribution function is often sufficient; this observation calls for the elaboration 
of a theory establishing the evolution of the reduced distribution, bypassing the evaluation of the full 
probability measure. 

A formal equation for the evolution of the mean occupation numbers Ni{v, k) is obtained by taking the 
expectation of the microdynamical equation ri{-, A; + 1) = £ri{-, k): 



N(-,fc + l) = E £ri{-,k) 

For the _Pi?— dynamics this equation reduces to 

Ni{Y + Ci,k + l) = E 



A; = 0, 1, 



(52) 



ttGII j 
ttGII j 



(53) 
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where we have used the microdynamical equation (27), the independence of the random vectors ^ and 
the Boolean field ri(-,k), and E[^^] = p^. Expression(53) provides a set of closed equations governing 
the evolution of the mean occupation numbers. These equations exhibit similarities with the continuous 
Boltzmann equations of statistical mechanics: they have the same structure and are obtained in a sim- 
ilar way; they are therefore referred to as the lattice Boltzmann equations of the model. Note however 
that lattice Boltzmann equations (53) exhibit two particular features: they are linear and are obtained 
without any approximation. These features arise from the absence of real interactions between particles: 
propagation and mixing are state independent channel permutations. Indeed, the only "interactions" 
between particles are in the velocity randomization because the exclusion principle does not allow inde- 
pendent changes of velocities of different particles on the same node. The correlations so produced do 
not preclude the derivation of an exact closed equation for the evolution of the mean occupation numbers 
but such correlations manifest themselves as soon as the evolution of two - or many - body functions 
is considered. For such many-body functions, it is still possible to obtain an exact closed equation but 
this cannot be achieved by a simple factorization of the expectation values^. It should also be noticed 
that the Boltzmann equations (53) have been obtained from the microdynamical equations written in a 
form which makes it clear that the dynamics proceeds by state independent channel permutations (27). 
When the automaton rules do not reduce to state independent channel permutations, the microdynamical 
equations must be written as in (29) and approximations must be used to obtain a set of closed equations 
for the densities Ni(r, k), a situation which will be discussed in the context of reactive systems. 



H. Macroscopic behavior 

To simplify the discussion, we consider the 4-equi-rotation velocity randomization operator (17) which 
transforms independently and simultaneously each node by performing a random rotation of the ve- 
locity configuration. When this velocity transformation is used in the _Pi?— dynamics, the Boltzmann 
equations (53) take the explicit form: 

#i(ri + I,r2,t + 1) = ^ [Ni{ri,r2, k) + #2(^1, r2, k) + #3(^1, r2, k) + #4(^1, r2, fc)) , 

#2(ri, r2 -Fl, t -Fl) = i (#2(ri, r2, k) + #3(^1, r2, k) + #4(^1, r2, k) + #i(ri, r2, fc)) , 

Nsiri - I,r2,t-M) = ^(N3{ri,r2,k) + N^{ri,r2,k) + Ni{ri,r2,k) + N2{ri,r2,k)^ , 

#4(ri, r2 - 1, t -Fl) = i (#4(ri, r2, k) + #i(ri, r2, k) + #2(^1, r2, k) + #3(^1, r2, fc)) . 

(54) 

This set of linear, finite-difference equations is conveniently solved in Fourier space. Introducing a Fourier 

mode"'^'^ 

#i(r, k) = Ni exp(ik • r) A\ (55) 
into (54) we obtain a linear algebraic system for the #j 's: 

4 

J2MijNj=0, i=l,...,4, (56) 



^Indeed, the BBGKY hierarchy is completely degenerate because equations (26) and (27) are linear in the 
variables rj, (r, k). 

^"Fourier indices will always be written in bold face and their components will be designated by numerical or 
roman subscripts; so there should be no confusion with the discrete time symbol k. 
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/ 1/4- Aexp(iA;i) 1/4 



M 



v 



1/4 1/4 
1/4 1/4- Aexp(iA;2) 1/4 1/4 

1/4 1/4 1/4- Aexp(-iA;i) 1/4 

1/4 1/4 1/4 1/4- Aexp(-iA;2) 



(57) 



Non-trivial solutions follow from the condition det M = and we obtain a fourth order polynomial 
equation for the damping factor A: 



detM = A^(A-^^^^^^^il±^^) = 



(58) 



The solutions A'-J\\i)(j = 1, ...,4) and the 



correspou' 



ding vectors N(-')(k) are given by-'^"'^ 



A^^\k) = 1/2 ( cos(A;i) + cos(A;2) ) ; N(i)(k) 



A(2)(k) = ;N(2)(k) 



A(3)(k) = ;N(3)(k) 



A(4)(k) = ;N(4)(k) 



f exp(-iA;i) 
exp(-iA;2) 
exp(iA;i) 

V exp(iA;2) 

/ 1 



-1 

V 

/ 
1 



V-1 

-1 
1 

V-i 



(59) 



(60) 



(61) 



(62) 



Since these vectors are linearly independent for all values of k, it follows that a general solution of the 
Boltzmann equation (54) is given by 



4 Li L2 



Ni{r,k) = J2J2 J2 a^^'\k)#}^')(k) exp f 27ri--iri j exp f 27ri-4r2 J M(j')(k) 

j = l ki=0k2 = ^ ^ ^ ^ 2 / V 

for a Li X L2 (periodic) lattice, and by 

4 /• 7r t'TT 

Ni{r, k) = Y. '^^^ '^'^2 a(j')(k) n''/\\l) exp(ik • r) {A^^\\L)f , 

j — l-J-n J -n 



(63) 



(64) 



for an infinite lattice, where the a(-')(k) are constants fixed by the initial condition. 

The macroscopic behavior of the system is determined by the dynamics of the long wave length Fourier 
modes 



exp (ik • r)(A(j')(k))*^N(j')(k) , |k| ^ 



(65) 



i.e., by the asymptotic properties of the damping factors j4'--'-'(k) when |k| 0. From the expressions (59)- 
(62) we have 



^^Here the kernel corresponding to A(k) = is spanned by the vectors (49) and (50), an appropriate choice as 
the components in this basis have a physical interpretation . 
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AW(k) = l-i|k|2 + 0(k4) , (66) 

and 

A(2)(k) = A(3)(k) = A(4)(k) = . (67) 
Two different types of temporal beliavior appear: 

1. Tlie damping factor A'--^'^ converges to 1 when |k| goes to zero ; i.e., the corresponding solutions 
relax arbitrarily slowly when the wavelength becomes very large. In the limit |k| = 0, A'--^'^ = 1 and 
the solution is time independent. 

2. The other damping factorsA^-* ) (j = 2,3,4) are equal to zero and the corresponding solutions relax 
to zero in one time step independently of the value of |k|. 

Consequently, when |k| is sufficiently small, there is one macroscopic mode persisting over long times (A ~ 
1), and three microscopic modes (A = 0) which are unobservable on macroscopic time scales. 

When |k| 0, the damping factor j4'-"'^-'(k) corresponding to the persisting mode is close to one and 
can be expressed as the exponential of an equivalent continuous damping rate w(k) 

A(i)(k) = exp(w(k)) , (68) 

or 

w(k) =ln(A(i)(k)) , (69) 



which, with (66), yields 



c.(k) = -i|kp + 0(|kn . (70) 



Equation (70) is the dispersion relation of the diffusion equation 
with 



o'^ 

p(^,t) = D—p(^,t) , (71) 



D = \. (72) 

In this sense, the long wavelength persisting modes in the LGA are diffusive. 

Notice that an isotropic dispersion relation is obtained to order 0(|kp) where the orientational dis- 
creteness of the lattice is not visible. However, anisotropic terms appear at order Odkj'*) where the fine 
structure of the lattice emerges. This result is not specific to the particular model considered and holds 
for a large variety of systems when discrete isotropy is present at the microscopic scale. 

A physical interpretation of the diffusive modes is obtained from the component of N(^)(k) in the basis 
{Nm, Nji, Nj2, Ng} (48)-(50). Expanding these components to first significant order in k yields 

4 + 0(|kn, (73) 
-2ifci+0(|k|3) , (74) 
-2ifc2 + 0(|k|3) , (75) 
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-2iki + 2ik2 + 0{\kf) , 



(76) 



respectively. We observe that the projection on the mass vector is dominant when |k| 0. Hence, the 
long wavelength persisting mode of the _Pi?— dynamics describes the diffusion of the mass density. 

It is worth mentioning that the i?_P— dynamics produces slightly different modes because, in that case, 
the state of the system is always observed after the velocity randomization. Indeed, when the velocity 
randomization is the 4-equi-rotation operator (17), the mass current density and the g— density are 
projected to zero before they are observed. On the contrary, in _Pi?— dynamics the state of the system is 
always observed after propagation, a step during which the channel densities propagate along the lattice 
links and produce local (/—densities and mass current densities where density gradients are present. 

I. Conserved quantities and spurious modes 

In the previous subsection, the time evolution of the long wave length modes was determined. Here we 
obtain a complementary characterization of the dynamics by considering the undamped modes and their 
spatial properties. 

The undamped modes of the _Pi?— dynamics with the 4-equi-rotation velocity randomization (17) are 
obtained from (59)-(62) by solving |j4(-')(k)| = 1 for j and k. Two solutions are found"'^^: 

1. For j = I and k = we have j4'-"'^-'(k) = 1 and the corresponding eigenvector N'-"'^-'(k) is the mass 
vector Nm (see (59)). We recover mass conservation. 

2. For j = 1 and fci = A;2 = tt we have j4("'^)(k) = —1 and the corresponding eigenvector N("'^)(k) is 
along the mass vector Nm (see (59)). We have an inhomogeneous mode oscillating with period 2. 
The corresponding invariant is the difference of mass on odd and even nodes (the parity of a node r 
is defined as the parity of the sum of its coordinates ri -|- r2). It is a cyclic invariant which changes 
sign at each time step. In the lattice gas automata literature, this is known as the checkerboard 
parity invariant. Invariants which, like the checkerboard parity invariant, do not have a counterpart 
in real systems are called spurious invariants, and the corresponding slowly decaying modes are 
called spurious modes [20], [21]. 

It should be noted that the wave vector k = (tt, tt) is not compatible with the boundary conditions 
when at least one of the system lengths Li or L2 is odd. In that case the checkerboard parity 
invariant is absent in the dynamics but very slow decays are observed for modes corresponding to 
wave vectors close to k = (tt, tt) ; their decay time is of the order of the square of the system size. 

A geometric interpretation of the checkerboard parity invariance can be given. If in a system with 
Li and L2 even, the lattice nodes are painted as a periodic checkerboard, each particle has a trajectory 
visiting alternatively a node of each color. As a result two particles on differently colored nodes will never 
interact and the cellular automaton universe consists of two totally independent subsystems corresponding 
to alternate colors at alternate time steps. 

If Li or L2 is odd, the lattice cannot be painted globally as a periodic checkerboard, but this can 
still be done locally. As a result two particles occupying neighboring nodes of different colors cannot 
interact before their mutual distance reaches at least one half of the smallest odd size of the system. This 
explains why the checkerboard parity mode becomes a very slow mode, with a decay time depending on 
the system size, when Li and/or L2 is odd. 

The cherkerboard parity invariance involves a microscopic wave length and therefore its infiuence 
on the macroscopic behavior can be thought to be negligible. This is indeed true when fiuctuations are 
ignored and when the Boltzmann equations are linear and obtained without any assumption. In this case, 
the modes around k = (tt, tt) are weakly excited by macroscopic initial conditions and the Boltzmann 



^^The same solutions are found when dynamics is considered. 
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equations guarantee that they remain so at any later time. In addition there is no coupling between 
modes, which protects the diffusive modes from the influence of the other modes. 

When reactive processes are included in the dynamics, it is not possible to derive a Boltzmann equation 
in a rigorous way because reactive interactions couple the Ni dynamics to correlation functions of arbi- 
trarily high order; in addition the Boltzmann equations are then usually non-linear. For these reasons, 
a mode initially excited with a very small amplitude is not guaranteed to remain so when the system 
evolves. Indeed, a weakly excited mode can be amplified to the macrosopic level when spontaneous 
symmetry breaking occurs. In this way, spurious modes can emerge at the macroscopic level even when 
they are absent from the initial condition; once excited, they can remain in the system for a very long 
or infinite time. Consider for instance a system with the checkerboard parity invariance, and suppose 
that this system is initially prepared in a microscopic state compatible with a smoothly varying density 
field Ni(r). If a spontaneous symmetry breaking occurs, it usually arises differently in the two checker- 
board subsystems because they are completely disconnected. As a result, the density field Ni(r,k) can 
evolve very differently in the two subsystems which makes the amplitude of the checkerboard parity mode 
grow from a microscopic fiuctuation to the macroscopic level. If the two subsystems are not considered 
separately, unphysical results can be obtained. If one wants to keep the two checkerboard subsystems 
unseparated, one can remove the checkerboard parity invariance by introducing at each node additional 
"channels" for rest particles (i.e. particles which do not propagate but are included in the velocity - or 
channel occupation - randomization). Indeed, the rest particles couple the checkerboard subsystems with 
one another, and as a result, the undamped checkerboard mode gives way to a decaying mode with a 
decay time depending on the number of rest particles per node and the efficiency of the mixing between 
rest and moving particles; with a small number of rest particles per node this time can be set to the order 
of an elementary time step. 

Since spurious modes can contaminate the physically interesting macroscopic behavior, they must be 
identified and, as far as possible, eliminated. If some spurious invariants remain in the dynamics, it 
is important to check that they do not affect the macroscopic behavior. It should be emphasized that 
spurious invariants are frequent in lattice gas automata because of the discrete nature of the models. As 
an example, the 2-eqm-rotatton velocity randomization involving only two different rotation angles 



that was used in the early developments of reactive lattice gas automata, yields up to seven spurious 
invariants when it is combined with propagation in a RP— or _Pi?— dynamics; a detailed description is 
presented in [22]. 

Other schemes for removing the checkerboard invariant that do not involve rest particles may also be 
constructed. For instance, the propagation and velocity randomization operations can be modified to 
couple the two sublattices or the nature of the lattice can be changed. Some of the simulations reported 
in the applications (Sec. VI to IX) were carried out on hexagonal lattices to avoid such spurious invariants. 



The long time behavior of the channel densities Ni(r,k) obtained in the previous section is an im- 
portant result but an incomplete characterization of the asymptotic dynamics since it does not provide 
any information about correlations. Here we show that a complete characterization of the asymptotic 
statistical state of the system can be obtained when the dynamics is based on the 4-equi-rotation velocity 
randomization operation and when the system has finite size. Only the main lines of reasoning will be 
sketched and the discussion will be restricted to a finite rectangular system with even dimensions Li and 
L2, and with periodic boundary conditions; a detailed discussion is presented in [22]. The important 
hypothesis is that the system is finite. One can proceed along the same lines for other geometries, but 
it should be emphasized that different geometries can produce a different asymptotic behavior because 
the presence of a global checkerboard parity invariant is very sensitive to the boundary conditions. Here 
Li and L2 are both even and therefore the system has the checkerboard parity invariance. The analysis 
given in [22] is as follows: 




with probability 1/2 , 
with probability 1/2 , 



(77) 



J. Equilibrium states 
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1. One first considers the evolution of the automaton at even times 

ry(r,fc) :fc = 0,2,4,6,...| ; (78) 

this process is a stationary, finite, Markov chain. 

2. Noting that the total number of particles remains constant separately in the two checkerboard 
subsystems, one splits the initial chain into reduced chains corresponding to fixed numbers of 
particles on odd and even nodes (at even times the two checkerboard subsystems can be identified 
with odd and even nodes respectively). 

3. One shows that the reduced chains so-obtained are aperiodic; this result is obtained by showing 
that the random sequence P o R o P o R reduces to the identity with non-zero probability. 

4. One also shows that the reduced chains are irreducible: starting from a given lattice configuration 
with Me and Mg particles on even and odd nodes, respectively, the dynamics can lead to any other 
lattice configuration with the same numbers of particles in the two subsystems. 

5. The above properties allow one to make use of a classical theorem (see for instance [23]) stating 
that a Markov chain which is finite, aperiodic and irreducible 

(a) has a unique invariant measure ; 

(b) is ergodic ; 

(c) is mixing (in the sense that any initial probability distribution converges towards the invariant 
measure). 

6. One shows that the probability measure which assigns the same weight to any lattice configuration 
with Me and Mg particles on even and odd nodes, respectively, is an invariant measure of the 
corresponding reduced chain. From the abovementioned theorem, this probability measure is the 
unique invariant measure towards which all initial measures converge. 

7. Having obtained the invariant measure at even times one determines what this measure becomes 
at odd times. 

This yields the following result. Suppose we have a dynamics based on the 4-equi-rotation velocity 
randomization combined with the propagation in a RP— or _Pi?— dynamics. Suppose also that the system 
resides on a finite rectangle of L1L2 nodes with Li and L2 even and with periodic boundary conditions. 
Suppose finally that the initial configuration ?y(-,0) has Me and Mg particles on even and odd nodes 
respectively. Then, the probability V(ri(-,k) = s(-)) to find the system in state s(-) at time k obeys the 
following relations 

"-'""'"^S^^S^j^-^^^-""^' ifstates(.)hasM.and 

Mo particles on even 



lim r{v{-,'2k) = s(-)) 

k^co 



and odd nodes, respec- 
tively ; 



(79) 



lim r{v{-,'2k + l) = s(-)) 

k^co 



, otherwise. 



(LiL2/2)!(LiL2/2)! ^^^^^ ^\ > ^"e 

and Mo particles on 
odd and even nodes, 
respectively ; 

, otherwise. 



(80) 



This provides a full characterization of the asymptotic behavior of any initial probability measure. 
Three consequences of these results are particularly important: 
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1. When a RP— or _Pi?— dynamics is based on the equi-rotation velocity randomization, the dynamics 
has exactly two independent global invariants: the total mass in each checkerboard subsystem. 

2. In each checkerboard subsystem the equilibrium state is fully characterized by a single parameter: 
the density per channel. 

3. In a large system at equilibrium, the occupation Boolean variables rii(r) of a small number / of 
channels (/ ^ L1L2) can be considered as almost mutually independent. 

K. Observation of diffusion 

In the previous subsections we showed how a probabilistic description of the automaton leads to 
diffusive behavior: the long wavelength modes which survive over long times are characterized by a 
dispersion relation whose leading term is diffusive: w(k) cx |kp. Here we briefly discuss the conditions 
that must be fulfilled to observe diffusive behavior in an automaton simulation. 

A first requirement is that the simulation must be performed on a large length scale so that the 
dispersion relation of the excited modes is close to w(k) cx |kp. Indeed the dispersion relation is never 
exactly diffusive (see (70) ) and the length scale depends on the required level of accuracy. Since the 
precision is frequently limited by other factors, wavelengths of the order of 20 to 80 lattice links are 
in practice often sufficient. Of course, long length scales imply long time scales for observing diffusive 
behavior; this is because, at a constant value of D, time scales like the square of the length. 

Since it is the density that diffuses (and not the particles which only perform random walks), diffusive 
behavior is seen only when the densities Ni{v, k) can be inferred from the random variables ?yi(r, k). The 
most direct way to do so is to consider an ensemble of rir independent realizations of the dynamics 

{ryf)(r,fc) : = 1, 2, . . . , , (81) 
and estimate the probability Ni(r, k) by the mean 

-Y^rft\v,k). (82) 

h = l 

Unfortunately, this method does not generalize nicely to reactive systems where important phenomena — 
e.g. spontaneous symmetry breaking — occur in single realizations of the dynamics. In these situations, 
performing an average over an ensemble of systems is still possible but not very relevant. For those systems 
the density p in the phenomenological reaction-diffusion equation and the lattice Boolean variables rii{v, k) 
must be connected in single realizations of the dynamics. This is accomplished by spatial and/or temporal 
averaging of the Boolean variables rii{v, k). However larger systems may be necessary to obtain a correct 
spatial resolution. 

L. Evaluation of the diffusion coefficient 

When performing simulations, it is quite useful to have a general expression for the diffusion coefficient 
in terms of the probabilities governing the velocity randomization. Then the diffusion coefficient can be 
given particular values by tuning the probabilities. For general velocity randomization and/or complex 
lattices, it may be difficult to obtain an explicit expression for D from the Boltzmann equations since 
this requires the solution of a generalized eigenvalue problem. 



^■^In addition, a spatio-temporal averaging acts like a filter which modifies the observed diffusion coefficient. This 
effect can be important when the wavelength is not very large compared to the typical length of averaging. 
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Here we derive a general expression for the diffusion coefficient by evaluating the time dependence of 
the mean square displacement of a tagged particle and by using Einstein's formula: 

where d is the dimension of space and (T^(k) is the average squared displacement performed by a tagged 
particle during a time interval A;. In a lattice gas automaton one must specify how a particle is tagged 
in order to use (83). Indeed, different tagging rules are possible and therefore different values of D can 
be obtained with Einstein's formula. This stems from the fact that the diffusion coefficient associated 
with a single tagged particle does not necessarily coincide with the diffusion coefficient that is relevant 
for mass transport in a many-particle system, which is what we are interested in. Therefore we should 
find a tagging rule such that both diffusion coefficients have the same value. However for the automaton 
with state independent velocity randomization, the diffusion coefficient is independent of the density and 
consequently we may as well consider a system with one single particle thus bypassing the tagging rule. 

In order to illustrate the method, consider the velocity randomization defined by (17) but where the 
four possible rotations can have different probabilities 

{Sj (r, k) with probability po 

Si+i(r,k) with probability pi 

Sj_l_2(i', fc) with probability p2 

Sj_l_3(r, k) with probability pi 

with Po + 2pi + P2 = 1- Note that the probabilities for the rotations 7r/2 and 37r/2 are equal according 
to our general assumption that the velocity randomization is invariant under the node symmetry group. 
Let 

{ri,r2,r3, . . .,rfc, . . .} (85) 

rk e {ci, i = 1...4} . (86) 

be the succesive displacements of the tagged particle. This sequence is a stationary Markov chain whose 
transition probability measure 

Pi J =V(rk+i = Cj\rk = Ci); i,j = l..A, (87) 

is given by 



/ Po Pi P2 Pi 

Pi Po Pi P2 

P2 Pi Po Pi 

V Pi P2 Pi Po 



(88) 



Except for special cases, this Markov chain is irreducible and aperiodic so that it has a single invariant 
measure towards which any initial measure converges; this is the uniform distribution. When the chain is 
reducible or aperiodic, a velocity related spurious invariant exists and the dynamics should be modified 
to remove it. Let 



d(fc) = ^r„, (89) 



be the total displacement after k time steps and assume that the initial condition is the invariant measure 

P(ri = Ci) = l/4, Vi; (90) 

by definition, the same distribution applies at any later time. Under these conditions, the expectation of 
d(k) is equal to at any time and 



25 



,[d{k)] = m\k)] = k + 2j2 E[rn' • r„] , (91) 



n' <^n 



which reduces further (c/, initial conditions and stationarity) to 

k-l 4 4 

2 



var[d(fc)] = + i ^(fc - n) c, . c, [P,,T . (92) 



n = l z = l j=l 

This expression can be rewritten as 

k-l 

2 



^ k-l 

[d{k)] = k+-Y,{k-n)Tr{GV-) , (93) 



n = l 



where the matrix G is Gij = Cj- • c j . The trace Tr(GP*^) on the r.h.s. of (93) is easily obtained in a 
basis where P and G are simultaneously diagonal (see Table I). The sum in (93) can then be evaluated 
and Einstein's formula (83) yields 

J=l ^°-^^ + ^ . (94) 
4 P2 - Po + 1 

For additional information on diffusion in lattice gas automata see Ref. [24]. 



Vectors 


Eigenvalues of P 


Eigenvalues of G 


(1,1,1,1) 


1 





(1,-1,1,-1) 


Po + P2 - 2pi 





(1,1,-1,-1) 


PO - P2 


2 


(1,-1,-1,1) 


PO - P2 


2 



TABLE I. Eigenvalues and eigenvectors of G and P. 
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IV. REACTION-DIFFUSION 



In a reactive LGA, the dynamics is built from three basic operators: the velocity randomization operator 
R, the propagation operator P, and the chemical transformation operator C , as described in Sec. II. At 
first sight, the chemical transformation operator C can be viewed as a simple extension of the velocity 
randomization operator R: both operators modify each node configuration independently by a stochastic 
transformation specified by a transition probability matrix. The difference is that R conserves the total 
number of particles while C does not. However, a subtle feature makes these two operators very different: 
when the velocity randomization is applied to a node configuration, the new average populations Ni can 
be predicted solely on the basis of the previous average populations (see (53)). This is not the case for the 
chemical transformation where the full joint probability distribution of the different channels is needed 
to predict how the average populations Ni are transformed. The difference arises because the occurrence 
of a particular reactive event a ^ /3 requires exactly a particles on the node before the reaction, a 
condition whose probability depends not only on the average population in each channel but also on 
correlations. As a result, a closed equation for the average populations Ni(r, k) cannot be obtained for a 
generic reactive LGA unless one considers particular regimes where approximations can be introduced to 
disconnect the dynamics of average populations from correlations. In this section, we consider the regime 
where reactions are infrequent and the Boltzmann approximation can be used to derive kinetic equations 
for the average population variables. The validity of these equations and how they lead to a macroscopic 
reaction-diffusion equation is discussed. 

A. From microdynamics to Boltzmann equations 

For simplicity, we restrict the presentation to the dynamics generated by an evolution operator of the 
form 

£ = PoCoR. (95) 

Before we begin the discussion of this dynamics, it is important to specify how a change in particle 
number a ^ /3, selected according to the transition matrix _P(a,/3), is translated into a configuration 
change s ^ s'. Each time particles are created or annihilated by the chemical transformation operator 
C , a modification is also produced in the velocity distribution. Unless models are designed to account 
for specific changes of the velocity distribution by reactive collisions, the simplest ansatz is to completely 
randomize the velocity distribution and neglect the reactive contributions to the diffusion coefficient 
However, this ansatz makes unnecessary changes in the velocity distribution and, therefore, does not 
minimize the reactive contributions to the diffusion coefficient. These contributions can be reduced by 
performing the transitions a ^ /3 in the following way: 

• When a non-reactive event a ^ a is selected, the initial configuration is not modified: s ^ s. 

• When a creation event is selected a ^ /3 (a < /3), the particles already present before the reaction 
are unaffected while new particles are created on empty channels (all combinations being equally 
probable). 

• When an annihilation a ^ /3 (a > /3) is selected, empty channels remain unaffected and particles 
are removed from initially occupied channels (all combinations being equally probable). 



^*For the particular models discussed in this section, the chemical transformation does not modify the diffusion 
coefficient because C is always combined with a velocity randomization R which reshuffies completely the velocity 
distribution. Models where the velocity randomization is not very efficient (R(s, s') with high probabilities on the 
diagonal) should be treated with caution. 
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These rules determine completely the matrix elements C(s, s') in terms of the particle number transition 
probabilities P(a,l3). The result is shown in Table II where it can be checked that 

4 

P{a,l3)= ^ C(s,s'), V s such that = a . (96) 

s'es «=i 

Equation (96) expresses the fact that any configuration s with a particles has a total probability P(a, [3) 
to be transformed into a configuration s' with [3 particles. 





0000 


1000 0100 


0010 0001 


1100 1010 


1001 0101 0011 0110 


1110 


1101 1011 0111 


1111 


0000 


-Poo 


Poi/4 Poi/4 Poi/4 Poi/4 


-P02/6 P02/6 P02/6 P02/6 P02/6 P02/6 


P03/4 P03/4 P03/4 P03/4 


-P04 


1000 


Pio 


Pii 





P12/3 P12/3 P12/3 


P13/3 P13/3 P13/3 


Pii 


0100 


Pio 


Pii 





P12/3 


P12/3 P12/3 


P13/3 P13/3 P13/3 


Pii 


0010 


Pio 





Pii 


P12/3 


P12/3 P12/3 


P13/3 


P13/3 P13/3 


Pii 


0001 


Pio 





Pii 





P12/3 P12/3 P12/3 





P13/3 P13/3 P13/3 


Pii 


1100 


P20 


P21/2 P21/2 





P22 





-P23/2 P23/2 


P2i 


1010 


P20 


P21/2 


P21/2 


P22 





P23/2 


P23/2 


P2i 


1001 


P20 


P21/2 


P21/2 





P22 





P23/2 P23/2 


P2i 


0101 


P20 


P21/2 


P21/2 





P22 





P23/2 P23/2 


P2i 


0011 


P20 





P21/2 P21/2 





P22 





P23/2 P23/2 


P2i 


0110 


P20 


P21/2 P21/2 





P22 


P23/2 


P23/2 


P2i 


1110 


P30 


P31/3 P31/3 P31/3 


P32/3 P32/3 


P32/3 


P33 





P3i 


1101 


P30 


P31/3 P31/3 


P31/3 


P32/3 


-P32/3 P32/3 





P33 


P3i 


1011 


P30 


P31/3 


P31/3 P31/3 


P32/3 P32/3 P32/3 





P33 


P3i 


0111 


P30 


P31/3 P31/3 P31/3 





P32/3 P32/3 P32/3 





P33 


P3i 


1111 


PiO 


P41/4 P41/4 P41/4 P41/4 


P42/6 P42/6 P42/6 P42/6 P42/6 P42/6 


P43/4 P43/4 P43/4 P43/4 


Pii 



TABLE II. Expression for the matrix elements C(s, s') in terms of the elements of P(a,fi) which are written 
as Paf) in the Table entries. The first column and row indicate the relevant configurations before and after the 
chemical transformation, respectively s and s', represented as 4-bit words (i.e. s = (si S2 S3 Si)). Configurations 
with same numbers of particles are grouped together, and the groups are separated by vertical and horizontal 
lines. 
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In the evolution operator (95), the velocity randomization and the chemical transformation are applied 
in sequence and it is convenient to define a collision operator which combines the two operations: 

n = RoC, (97) 
= CoR. (98) 

Both R and C modify the lattice configuration by independent node transformations, and this property 
extends to TZ which can be described by a local transition probability matrix with elements 

{7^(s,s') : s,s' e 5} , (99) 

giving for each pair of configurations s and s' the probability TZ{s,s') that s is mapped onto s' by TZ. 
From the definitions (97) and (98) this transition matrix is simply the matrix product of the transition 
matrices of C and R: 

s"es 

= ^ C7(s,s")i?(s",s') . (100) 
s"es 

To derive the microdynamical equations corresponding to the dynamics generated by (95), we write 
the evolution operator as 

£ = Pon, (101) 

and we notice that this is identical to £ = P o R except that the transition probability matrix of R must 
now be replaced by TZ(s, s'). We can therefore use the same form of the microdynamical equations as for 
the i?_P— dynamics in Sec. Ill (see (29)): 

4 

ry,(r + c„k+l) = ^6s'(r, k) s', ^^iv, k){l - ry,(r, fc))(i-^^) , (102) 

s,s' j = l 

where we change the probability distribution of the Boolean matrices [^ss'(^) ^)]- 

1. for each configuration s there is one and only one configuration s' such that 
^ss'(r,fc) = 1, 

2. the probability of the event ^ss'(^) A;) = 1 is given by TZ{s, s') , 

3. random matrices [^ss'(^) ^)] different nodes and/or different times are mutually independent. 

The interpretation of the random matrices remains the same: a node configuration ri{v, k) is transformed 
by TZ into the only configuration s' such that ^-qs'i^, k) = I. 
Consider the identity 

4 

^^ = j2s,i[^;^{i-^jf-^^^ , (103) 

s j = l 

which follows from the fact that 

4 

n^yf (r,fc)(l-ry,(r,fc))(i-^^) , (104) 
i=i 

plays the role of a "matching indicator" : its value is equal to one if the configuration at node r at time 
k matches the configuration s (i.e. if T7(r, A;) = s) and zero otherwise. Equation (103) together with 
J2s' iss' = 1 yields 
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Combining (102) with (105), we obtain 

4 

ry,(r + c„ k + 1) - ry,(r, fc) = " 6s' (r, A;) 7j^^(r, k){l - ry,(r, k) f-'^'> 



(105) 



(106) 



which is equivalent to (102) and will prove to be a convenient form for later applications. 

A formal equation for the evolution of the mean occupation numbers Ni(r, k) is obtained by taking the 
expectation of the microdynamical equation (102) 



Ni{v + Ci,k+X) = E 



(107) 



Using the independence of the random matrices ^ and the Boolean field ri{-, k), and using E[^ss'(i', k)\ 
TZ(s, s'), we simplify this equation to obtain 



N,{r + c„k + l) = ^7^(s, s') s', E ^ m - ry,(r, kjf-^^^ 



This is not a closed equation for the mean occupation numbers Ni(r, k) since the expectation 



(108) 



n'?f(r,fc)(l-ry,(r,fc))(i-^^) 



(109) 



involves not only average occupation numbers but higher moments, refiecting correlations between dif- 
ferent channels at the same node. This is a consequence of the fact that the matching indicator (104) 
involves the occupations of all channels at node r at time k, so that the knowledge of the complete 
probability distribution of ri{v, k) is required to evaluate the expectation (109) for all the possible values 
ofs e S. 

At this point it is instructive to consider again the diffusive models presented in Sec. Ill in order 
to understand why the dynamics of the average populations is independent of the correlations in these 
models. There the decoupling arises because for each possible realization of the velocity randomization 
R, the post-randomization state of a channel w{v, i) at node r depends on the state of a single channel 
w{v,j) at the same node. The particular channel w{v,j) which determines how the channel w{v,i) is 
transformed depends, of course, on the realization of R but this is unimportant: the crucial point is 
that, no matter the result of the random selection, for each realization of R there is only one channel 
w{v,j) whose state infiuences the transformation of the channel w{v, i). In a reactive model this cannot 
be expected because the way in which a node is transformed by the chemical operator depends precisely 
on the number of particles at the node. 

In general, the expectation (109) cannot be expressed in terms of the average ocupation numbers 
Ni(r,k). However an important exception arises when the variables rii{v,k) are decorrelated: in this 
case (109) can be replaced by 



HA^f (r,fc)(l-^^,(r,fc))(l- 



(110) 



^Notice that in (109), the power in any of the variables r]t is equal to 1 or 0. 
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Therefore, if we assume that the variables rii{v,k) are decorrelated before each application of 7?. - this 
is the Boltzmann ansaiz - we can substitute (110) for (109) in the formal equation (108) and obtain a 
closed equation for the average occupation numbers: 

4 

N,{r + c„k+l) = ^7^(s, s')s', n - ^))^'"'^'^ > (111) 

or equivalently from (106) 

4 

Ni{r + ci,k + l)- Ni{r, k) = ^7^(s, s'){s'i - s,-) J] Np{r, k){l - Nj{r, kjf-'^^ . (112) 

s,s' j = l 

This set of nonlinear, finite-difference equations is analogous to the continuous Boltzmann equation of 
statistical mechanics; ^'^ for this reason (112) (or equivalently (111)) is referred to as the latttce- Boltzmann 
equation of the lattice gas. 

It is usually argued that the Boltzmann approximation should hold in the limit of small gradients and 
infrequent reactions 

C(s,s') < 1 for s 7^ s' , (113) 

since then reactive processes occur on nodes which support a near diffusive equilibrium distribution where 
the velocity channels are uncorrelated. In such a circumstance the rate at which the diffusive equilibrium 
is perturbed by the reactive process is much slower than the rate at which the system returns to a 
local diffusive equilibrium. However, the argument should be considered with caution because a typical 
relaxation rate is meaningless for diffusion in an infinite system. Indeed, the diffusion rate 

\uj\ = D\kW (114) 

depends on the length scale (|k|~"'^), and on large enough length scales, the diffusive relaxation is always 
slower than any reactive process. Therefore large scale phenomena can compromise the validity of the 
Boltzmann approximation even for vanishing reactive rates. These considerations imply that the 
reaction-diffusion equations follow only in a suitable weak coupling limit that considers both spatial 
gradients and the ratio of reactive to non-reactive relaxation times to be small. [26] In addition, it has 
been shown that the correlations which are omitted in the Boltzmann-level description are responsible 
for a number of physically relevant phenomena. For instance, the diffusion limit, where the rate constant 
takes the Smoluchowski form and is no longer kinetically dominated, has its microscopic origin in infinite 
sequences of ring collision events which are not considered in the Boltzmann equation. [27] 

Correlation effects due to incomplete mixing are especially evident in low dimensional systems and 
numerical simulations based on random walk models [28] have provided the basis for a large number of 
investigations addressing the question of the validity of the local diffusive equilibrium assumption with 
the aim of increasing our understanding of the new types of kinetic behavior that emerge when this 
hypothesis is not applicable [10]. 

The Boltzmann equation may also be inadequate when the microscopic randomness does not disappear 
at the macroscopic scale. Typical examples are spontaneous symmetry breaking or nucleation processes; in 
both cases, microscopic fiuctuations trigger persistent large scale phenomena. In spite of these limitations, 
the Boltzmann equations provide a valid description for many situations of practical interest. Examples 
where the Boltzmann equations apply and where they do not are discussed in the applications sections. 
Non-Boltzmann phenomena in lattice gas automata have been investigated in Ref. [29], [30]. 



^^This is most easily seen by noting that in (112) the r.h.s. is a coUisional term and by taking the continuous 
limit of the L h. s.. 

^'^The idea that the local equilibrium hypothesis could be inappropriate for some reactive systems was first 
suggested around 1950 in the context of exothermic reactions [25]. 
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B. From Boltzmann equations to reaction-diffusion behavior 



An alternative form of the Boltzmann equations, (111) and (112), is obtained when the average popu- 
lation variables Ni are expressed in terms of the mass density p(r, k), the mass current density j(r, k), and 
the q-density pq{v,k) (see (42), (43) and (51)). This transformation is interesting because it separates 
variables into slow and fast types. Indeed, the q-density and the mass current are not conserved in the 
velocity randomization operation and therefore these variables are expected to exhibit fast relaxation (in 
a few time steps). On the other hand, creation and annihilation of mass occur only via reactive processes 
which are assumed to be infrequent. When gradients are weak the relaxation of the mass density is 
dominated by reactive processes and takes place on a long time scale. If microscopic time scales are 
ignored, the fast variables can be eliminated to obtain a description which focuses on the slow dynamics 
of the mass density. How this elimination leads to a reaction-diffusion equation will not be discussed in 
full generality; instead, we consider the situation where the Boltzmann equations are linearized around a 
homogeneous steady state. 

1. Homogeneous solutions of the Boltzmann equations 

We show that the determination of a homogeneous and isotropic solution 

N^{v,k) = p,lA = c,, Vi= l,...4;fc = 0,1,2... ;re£ , (115) 

for the Boltzmann equation (111) is a problem equivalent to solving the algebraic equation 

/(c.) = 0, (116) 

where / is the macroscopic reactive rate (7) and Cg denotes the stationary value of the density of species 
X per channel, cx = Cg. Substituting Cg for Ni{v,k) in the Boltzmann equation (111) yields 

4 

c, = ^s:.7^(s,s')^c:^■(l-c.)^'-^^■^ • (117) 

Since C and R are invariant under rotations and reflections of the channels, the same invariance applies 
to TZ] ^® therefore the r.h.s. of (117) is independent of the index i: 

s,s' i = l j = l 

Expressing 7^(s, s') in terms of C(s, s') and i?(s, s') as in (100) we obtain 

S,S" j = l s'eS 8 = 1 

which reduces to 



^^By construction, the collision operator conserves the invariance property of R and C: 

Tl{s,s') = Tl{gs',gs) , Vs, s' G 5 , (118) 

for all node transformation g corresponding to a rotation or a reflection of the channels in real space. 
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= i E E ^(^' ri <'(^ - ' (121) 

s,s' z' = l i = l 

where we have used the conservation of mass during the velocity randomization: 

j2ij2^-)Ri^",^') = J2'''' ■ (122) 

s' i = l i = l 

Using the identity 

4 

= ^C7(s,s') (1 - N/f-'^^ , (123) 

s,s' i = l 

which follows from C(s, s') = 1, and summing the result over the index i, we can rewrite (121) as 

= I E E(^'- - ^') ^(^' ^') ri ^^^(1 - '^y'"''' ■ (124) 



Combining the different terms on the r.h.s. which correspond to pairs of configurations s and s' with a 
and l3 particles, respectively, (see (96) or Table II), we find 



= l/(^^) = E E l(/5 - ")^("> /?) (!) c?(l - cO'-" • (125) 



4 4 
a = 13 = 

This shows that the homogeneous isotropic solutions of the Boltzmann equations coincide with the fixed 
points of the phenomenological rate equation 

^ = i/(«), (126, 

written here in terms of the density per channel. 



2. Linearized Boltzmann equations 

If we assume that the system is close to a homogeneous isotropic steady state we can linearize the 
Boltzmann equation (112) about this state. Defining 

6Niir,k) = Niir,k)-c, , (127) 

where Cg is a solution of /(c^) = 0, we have to first order in 6Ni 

6Ni{r + Ci,k + 1)- 6Ni{r, k) = jCij6Nj{r, k) , (128) 

where 



= ^^,(N) 



(129) 

{Ni=c, :£=1,...,4} 



with 

4 

^,(N) = ^7^(s,s')(.:■-.,)^^f (l-^^-)^'"^^'^ • (130) 
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It is sometimes useful to write (128) as 

6Ni{r + ci,k + l) = Bij6Nj{r, k) , (131) 

with 

Bij = Sij + jCij . (132) 

The matrix [Bij] describes how small deviations from Cg evolve when the collision operator TZ is applied to 
a node. The main properties of the linearized collision operator [Cij] and of the matrix [Bij] are discussed 
in the next section. 



3. Properties of the linearized collision operator 
From (129) and (130) we have 

A, = E(s:--s,)^(s,s')(2s, -1) n C(l-c,y-'' 

which may also be written in the form 

A, = Es:•7^(s,s')s, n ^Hi-csf-'" 

-^.:.7^(s,s')(l-.,) n cJHi-c.)'-^" • 



(133) 



(134) 



The first term of the r.h.s. accounts for changes in Ni in a collision when there is a particle in channel 
w{v,j) before collision. The second term on the r.h.s. is the corresponding quantity when channel w{v,j) 
is empty before collision. 

As a consequence of the invariance of TZ under the group of rotations and refiections of a node on the 
lattice, the matrix [Cij] has necessarily the form 



/ Loo Loi Lo2 Loi 

Loi Loo Loi Lo2 

Lo2 Loi Loo Loi 

\ Loi Lo2 Loi Loo 



and in view of the definition (132) the matrix [Bij] has the same form 



(135) 



/ -Boo -Boi Bo2 Boi 

Boi Boo Boi Bo2 

Bo2 Boi Boo Boi 

\ Boi Bo2 Boi Boo 



(136) 



Now we show how the linearized collision operator [Cij] is connected to the macroscopic rate f(cx)- 
Summing the expression (133) over i and using (100) and (122) we obtain 



ea.=ee(^'--^'M«,«')(2s.-i) n ^j'd-c.)^ 

i = l s,s' i = l k,k^i 



(137) 



To arrive at (133) account is taken of the Boolean nature of Sj. 
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On the r.h.s. of this expression, we combine the terms corresponding to configurations s and s' with a 
and /3 particles respectively (see (96) or Table II). This leads to 



i — l 

In the sequel, we consider only those values of Cg which correspond to linearly stable fixed points of the 
rate equation (126), i.e. k is always positive. 

4. Solutions of the linearized Boltzmann equations 

In this subsection, we give the general solutions of the linearized Boltzmann equations for the dynamics 
based on the velocity randomization performed with the 4-equi-rotation rule defined in (17); in this case 
the matrix [Bij] has a form which allows analytical calculations. These manipulations parallel those given 
in (55) to (70) and are repeated here for the sake of clarity. When the velocity randomization (17) is used 
in TZ (see (97)), all the channels on a node are made statistically equivalent after collision. Therefore the 
matrix [Bij] takes the simple form 



f B B B B 

B B B B 

B B B B 

\B B B B 



(139) 



and the Boltzmann equations read 

4 

6Ni{r + Ci,k+l) = Bj2^Nj{r,k) , i=l,..A. (140) 

i=i 

Inserting in (140) the Fourier mode 

6Ni{r, k) = 8Ni exp(ik • y)A^ (141) 

we obtain a linear algebraic set for the 6Ni's 

4 

J2Mij6Nj=0, i = l,...,4, (142) 

i=i 

/ B -A exp(i ki) B B B 



M 



B B - A exp(i ^2) B B 

B B B- A exp(-i ki) B 

\ B B B B - A exp(-i ^2) 



(143) 



Non-trivial solutions follow from the condition rfetM = 0. Making explicit use of this condition, we 
obtain a fourth order polynomial equation for the damping coefficient A 



i3/i cos(A;i) + cos(A;2), 



Aei m = A-" {A - B ^) = . (144) 

The solutions A^j\\i){i = 1, ...,4) and the corresponding vectors (^N*--* -*(!€) are given by^*^ 



"'Again the kernel corresponding to j4(k) = is spanned by the vectors (49) and (50). 
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A'-^^k) = 25 ( cos(A;i) + cos(A;2) 



,5#(i)(k) 



(145) 



A(2)(k) = 



,5#(2)(k) 



/ 1 



-1 

V 



(146) 



A(3)(k) = 



,5#(3)(k) 



/ 
1 



V-1 



(147) 



A(4)(k) = 



,5#(4)(k) 



/ 1 

-1 

1 

V-1 



(148) 



Two different types of temporal behavior appear : 

1. The damping factors A'-J'^ (j = 2,3,4) are equal to zero and the corresponding solutions decay to 
zero in one time step independently of the value of |k|. This reflects the strong dissipation of the mass 
current density and the q-density as can be checked by projecting (146)-(148) onto (48)-(50). Notice 
that the "instantaneous" decay ^(j) = (j = 2,3,4) arises because the velocity randomization operators 
erase all information about the particle velocities in one single time step. 

2. The damping coefficient A^^) depends on k and its value is in general different from zero. 
Since 



we have from (132) and (138) 



45 = E^': 



iB = 1- K 



(149) 



(150) 



so that the damping coefficient reads 

AW(k) = (l-.)(^^^^5^^il±^^) . 

For small wave numbers and infrequent reactive collisions (i.e. k <^ 1), the damping coefficient j4(^)(k) 
is close to one and can be expressed as the exponential of an equivalent continuous damping rate w(k) 



(151) 



A(i)(k) = exp(w(k)) 



(152) 



;(k) = In (A(i)(k)^ 

= ln(l-«)- |k|V4 + 0(|k|4) 



(153) 
(154) 



Equation (154) shows that in this regime the dispersion relation is equivalent to that of the linearized 
reaction-diffusion equation 



d 

— 8cx{v,t) = -K6cxir,t) + DVHcxir,t) 



(155) 



with D = 1/4. In this equation, the field (5cx(r, t) = cx(r, t) — Cg can be identified with the fiuctuation of 
the mass density per channel on the lattice because (145) is oriented exactly along the mass vector (48). 
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V. REACTIVE LATTICE GAS 



The formalism required to properly describe a multi-species automaton dynamics in all generality 
and to carry out the statistical mechanical analysis which yields the reaction-diffusion equations is rather 
involved and will not be presented here. (The microdynamical formulation of the multi-species model can 
be found in [18].) However, it is easy to write a general algorithm suitable for simulation of the reactive 
dynamics without giving details of the microdynamical theory. This section provides the "recipe" for such 
an algorithm and can form the starting point for the development of computer codes for implementing 
the reactive lattice gas automaton. 

We consider reactive systems which are described phenomenologically by partial differential equations 
of the reaction-diffusion type 

^^^ = F{p{r,t)) + B-V'p{r,t), (156) 

where D is a diagonal matrix with positive diagonal elements. Here, each component pr (r = 1, 2, . . . , n) 
of the column vector p represents the local concentrations at space point r at time t of a corresponding 
species Xr involved in the reactive processes described by the reaction rate vector F(p). 

On a phenomenological level the reactive terms in (156) derive from the mass action rate law 

^ = F(p(t)) , (157) 

which follows from the reactive collision processes specified in a reaction mechanism. Therefore, F 
contains only polynomial contributions in the species concentrations. The second term on the r.h.s of 
(156) describes diffusive transport. 



A. Multiple species and exclusion principle 

Single-species models generalize easily to allow for different reactive species to move and react on a 
lattice. As previously, the positions of the particles are restricted to a regular lattice £, the time takes 
integer values, and the particle velocities take their values in a finite set of vectors compatible with the 
discretization of space and time 

There are many different ways to generalize the exclusion principle to account for the existence of 
different species. For instance one can allow different particles to be at a same node r provided they belong 
to different species and/or have different velocities. With this rule, each species is subject separately 
to its own single-species exclusion principle and a useful representation of the lattice is in terms of a 
"stack" of species lattices r = 1, • • where each species r resides on its own lattice where the 
single species exclusion rule applies. Thus, the lattice C may be decomposed into n species lattices 
C = Ci X ... X £„, with identical node labels r. This notional decomposition of C is useful both 
conceptually and in the mathematical formulation of the automaton. In this way non-reactive collisions 
can be carried out separately on each lattice Cj with different frequencies, mimicking the elastic collisions 
with the solvent, and reactive collisions couple the dynamics on all the species lattices. With this 
generalization of the exclusion principle, we shall see that the diffusion coefficients of different species can 
be tuned independently. 

It is straigthforward to extend the notion of channel. At each node r we define nm channels 
Wj{v,i), T = i = l,...,m where n and m are respectively the number of different species 



^^In some circumstances, it may be convenient to associate different sets of velocities to different species or 
restrict species to move on sublattices of C (some of wliicli can be identical and some can be C itself). We will 
not consider sucli generalizations liere. 
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and the number of different velocities in the model. A channel Wr(r, i) is said to be occupied if there is 
a particle Xr with velocity Cj- on node r (otherwise the channel is said to be unoccupied). 

A more restrictive exclusion principle is obtained when a channel cannot be occupied simultaneously by 
particles with the same velocity independently of their species (global exclusion principle). [31] It is also 
possible to apply the latter rule separately on different subsets of species: if two particles belong to the 
same subset the restrictive rule applies while if they belong to different subsets they are allowed to occupy 
the same channel. All the models so constructed are conceptually similar but differ mainly by (i) the 
class of reactive rates F(p(t)) to which they lead, and (ii) the flexibility they offer to tune independently 
the diffusion coefficients of different species. Here we restrict ourselves to the single species exclusion 
principle (In sec. X we consider some models without exclusion principle). 

B. Automaton rule 

The automaton rule is built from the composition of propagation, velocity randomization and chemical 
transformation operators: 

• propagation operator, _P^: Moves particles of species r a specified number of lattice units in the 
directions determined by their velocities to other nodes of the lattice. 

• velocity randomization operator, i?^: Randomizes the particle velocity configuration on the 
lattice Cr at each node independently of the others. The mixing operations, governed by the 
operators Rr, conserve the number of particles of species r; however, their momentum is not 
conserved locally. 

• chemical transformation operator, C: Is responsible for local chemical reactions among the 
species. Denoting by a = (aia2 • • - an) and /3 = {(iiji'z • • ■ Pn) the vectors specifying the particle 
occupancy at a node, where < aj,j3j < m for each r = 1, the reaction occurs with 
probability P{oL,(j), which depends only on the occupancy of the nodes r on £ and not on the 
velocity configurations. In order to fully specify the reactive transformation, the matrix P = 
[_P(a,/3)] must be completed by a rule detailing how particle number changes are implemented into 
channel occupation number changes. The precise knowledge of this implementation is important 
to properly formulate the automaton dynamics. However we will see that this level of description 
is not required to determine the phenomenological rate law F of an automaton dynamics; the 
knowledge of P is sufficient for that purpose. In this section we will therefore ignore how particle 
number changes are actually translated into node configuration changes. However we will make the 
assumption that the configuration of a species Xj remains unchanged in any reaction which does 
not change the number of particles of species Xj (a^ = j3j). In this way, we can neglect reactive 
contributions to the diffusion constants when reactions are unfrequent (see Sec. IV). 

Neglecting reactions for the moment, the (non-reactive) dynamics of the system is given by the compo- 
sition of free streaming and elastic collisions, PjoRj. To account for the possibly different elastic collision 
frequencies we may apply this composition of operators different numbers of times for each species in one 
total automaton time step. The general expression for the non-reactive automaton evolution is 

[State at time {k -|- 1)] = 5 [State at time k\ (158) 

with 

n 

£=]\{PrORrt , (159) 
r = l 

and Ij an integer. Combining this non-reactive dynamics with the chemical transformation C , we obtain 
an evolution rule for the full reactive dynamics: 

£={f{{PrORrt\c . (160) 
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C. Reaction probability matrix 



In order to complete the "recipe" for the rule construction some guidelines need to be given for the 
specification of the reaction probability matrix. We give a general strategy for the determination of 
P. In most circumstances the macroscopic equations of motion provide a good description of the basic 
phenomenology of non-equilibrium chemical rate processes. If the system is spatially homogeneous (well 
stirred) the appropriate equations of motion are the laws of mass action kinetics, while if the system 
is inhomogeneous diffusion terms can be used to augment the description leading to reaction-diffusion 
equations. Since the mass action rate law provides a mean field level of description of the local kinetics, 
it is natural to demand that the mean field approximation to the automaton dynamics correspond to the 
phenomonological rate law. We make use of this correspondence below. It is not difficult to derive an 
expression for the time rate of change of the average species concentrations from a mean field description 
of the automaton dynamics. We now summarize and formalize what was presented as a simple physical 
picture in Sec. II. 

Suppose that the system is spatially homogeneous and correlations are neglected. Then the probability 
to find a node in a configuration a is given by 



n 



^ ;c««(i-c.)™-«« , (161) 



since the particles are binomially distributed. In this equation = Pr/m, where m is the number of 
channels per node. Given that the probability of a transition from a configuration a to a configuration 
/3 is P{oL,(3) the rate of change of the average concentration of species r per node, pj, is 



p.(fc + l)-p.(fc)= ^(/3.-aOP(«,/3)n c:«(l-c.)™-«« (162) 

a K=i ^ ^ 

Here 

Pr{a) = Y,{pr-ar)P{a,l3) , (164) 
/3 

is the average value of the particle number change of species r when the input configuration is a. The 
rate law (163) describes the mean field reactive dynamics of the automaton particles. We now establish 
a correspondence between the fictitious automaton universe and the phenomenological description of 
reacting systems. 

Consider a general reaction mechanism consisting of r elementary steps 

v{Xi + + ■ --vlXn f v{Xi + 4X2 + ■ --^Xn, (165) 

where j = 1, • • • , r labels each step in the mechanism, which is completely characterized by the sets of 
stoichiometric coefficients = (i^^, 4^ ' ' ' j ^li) j — {4' 4'' ' ' ' ^n) ^^'^ ^^e rate coefficients {k±j : j = 
1, • • • , r}. The mass action rate law can be written in the form 



=1 




^ = }^(^^-^^)<;-fc, + Up?;. . (166) 



In order to relate the mean field automaton equation (163) and the mass action rate law (166) we assume 
that the concentration units in these two equations are the same and express the time in terms of the 
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automaton time unit, h. In these units the discrete-time version of the continuous-time mass action law 
(166) may be written as 

r n , n , ^ 

Prik + h) - prik) = hJ2 (4 -kj n P-^ + ^-i n f + ^e*') • (167) 

j=l I K=l K=l J 

We now determine the conditions on the reaction probability matrix that ensure the descriptions of 
the reactive dynamics given in (163) and (167) agree to order h. These equations can be made to have 
the same forms either by expressing each pr in (167) in the basis {c'^(l — c^)™"' : i = l,...,m} 
or, alternatively, expanding these factors in (163) to yield powers of Pr-^^ Choosing to expand in the 
{c\.(l — Cj)"^~^ : i = I, . . . , m} basis, we obtain 

r mm 

p,(k + h)-prik) = hj2ii.i-Di)j2---j2 

j=i ii=o ;„=o 




(168) 



Identification of the right hand sides of (163) and (168) yields a set of constraint conditions on the 
automaton reaction probabilities: 

r 

i=i 

n 

k., n 

K = l 

where (') = if a > h. These conditions on -Pr(«) ensure that the mean field description of the spatially 
homogeneous automaton dynamics is the same as the mass action rate law. Since Pt{ol) is an average 
value of the particle number change, (169) fixes only this particular moment; higher moments which 
refiect correlations of particle number changes are not determined by the mean field mass action rate law. 
There is no guarantee that the full automaton dynamics will yield such a mass action rate law. Indeed, 
one of the main interests in the automaton dynamics is that it contains correlations and fiuctuations that 
go beyond such mean field descriptions. Nevertheless, this scheme is very useful for tuning the reactive 
lattice gas parameters so that the automaton is likely to yield macroscopic behavior of a given type, an 
interesting feature as it is obtained from a mesoscopic description. 

Consider the case where the rate law is characterized by the set of rate constants {k±j : j = 1, • • • , r}. 
Since P{oL,(j) are transition probabilities, Pt{ol) as given by (164) is bounded: — < Pt{ol) < m — ctj. 
In view of the time scaling factor h in (169) the values of the Pj{ol) may be scaled to lie arbitrarily 
close to zero; however, their signs cannot be changed. Consider, for example, the extreme cases where 
= or m. Then, from (164), we have Pj{ol) > if = and Pj{ol) < if = w. The kinetic 
parameters on the right hand side of (169) must be such that these sign conditions are met. If the sign 



■^■^The number of channels per node m must be equal to, or larger than the highest power in pr in the mass-action 
rate law. 



-.;,{-*,n".'^(::a) (:)"'+ 
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conditions are met the scaling factor h can always be chosen such that P{oL,(j) is a probability. If the 
sign conditions are not met then the kinetic scheme characterized by a particular set of rate constants 
{^±j '■ i = 1) • • • ) J"} cannot be simulated directly by the automaton. 

Additional insight into the sign conditions on Pt{ol) can be obtained by examining the fluxes at the 
surfaces of the hypercube Q < pj < m within which the species densities must lie as a consequence of 
the exclusion principle. The fluxes on the boundaries of this cube must always be directed to its interior 
(or be tangent to the boundary). This follows from (163) and the sign conditions determined above for 
Pj{ol). Furthermore, one sees from the right hand side of (169) for = or m that one obtains the 
same sign conditions on Pj{ol). Thus, no additional constraints on the reaction probability matrix are 
implied by the flux condition. 

The elements of P are not uniquely determined by the set of macroscopic rate constants and the sign 
conditions and additional physical considerations that account for correlations in particle number changes 
can be used to construct a series of automaton models. This flexibility in the choice of P can be exploited 
to construct an automaton dynamics that goes beyond the mean held level. 

These formal considerations suggest a strategy for the construction of the reaction probability matrix. 

• Check that the set of rate constants in the mass action rate law is admissible by calculating the 
species fluxes at the boundaries. 

• Express the order h term on the right-hand side of (167) in the basis 

|n (rj^^^^l-^'^)™""" : «'^=0,...,m| (170) 

and identify each component with the corresponding Pj{ol). This yields (169). 

• For each value of a, flnd a set of positive values (possibly larger than 1) of _P(a,/3) satisfying (164). 
Physical considerations concerning correlations should be used to select the elements of P. 

• Considering all a, flnd the largest value of 

^P(«,/3). (171) 

Select a value for the scaling factor h in (169) such that the sum in (171) is less than unity. For 
each value of a the diagonal elements are then given by 

P(«,«)) = l- J2 (172) 

By construction P is a transition probability matrix. 

This completes the speciflcation of the algorithm for the construction of reactive lattice-gas automaton 
rules. There is a variety of ways in which the propagation and velocity randomization steps can be carried 
out and there is flexibility in the construction of the reaction probability matrix, as we have seen above. 



In such a circumstance reactive lattice gas models without exclusion can be used; these are briefly discussed 
in Sec. X. For some specific reactive schemes, the problem can also be solved by using the same lattice C for 
all species with an exclusion principle that forbids more that one particle to be in the same channel at the same 
time (instead of a stack of lattices Cr and an exclusion principle for each lattice). An important drawback of this 
latter method is the difficulty to assign different diffusion coefficients to different species [31]. 



41 



D. Applications 



In sections VI to IX we discuss applications to specific cliemical systems and examples of the construc- 
tion of rules will be given. 

We have argued that reactive lattice-gas automata can be used to investigate the dynamics of spatially- 
distributed chemically reacting systems at the mesoscopic level. The examples in the folowing sections 
have been chosen to illustrate a number of important spatial and temporal structures that are com- 
monly observed in reacting media. The emphasis in these examples is on the role of fluctuations on the 
far-from-equilibrium dynamics and on the schemes that are used to construct the automaton models. 
The phenomena that we consider include pattern formation and wave propagation in bistable chemical 
systems, spiral waves in excitable and oscillatory media, Turing pattern formation and the dynamics of 
chemical systems with deterministic chaos. For each case simulation results are shown and analysed. 

In a coarse-grained view of the dynamics the magnitude of the local particle number fluctuations 
depends on the size of the fluid volume element and the time scale on which the fluctuations are measured. 
The automaton reaction probability matrix P determines the rate at which chemical reactions take place 
at the nodes of the lattice and, therefore, controls the local particle number fluctuations that arise from 
such reactions. Thus the magnitude of the local particle number fluctuations and the correspondence 
between the automaton space and time scales and those of the real system are related topics which are 
intimately connected with the construction of the reaction probability matrix. In the course of discussing 
the applications we shall also describe how automaton models can be constructed that correspond to 
different space-time coarse grainings of the real reactive dynamics. Basically, if the automaton elementary 
cell volume represents a small fluid element and the automaton time step corresponds to a short time 
interval, the details of the reaction mechanism will be important and will determine the nature of the 
local fluctuations. However, if the automaton elementary cell represents a larger fluid volume and/or the 
automaton time step corresponds to a longer time interval, many reactions occur in a local cell in a given 
time interval and only the macroscopic kinetics embodied in the mass action rate law will be important; 
the details of the mechanism are then inessential, and indeed many reaction mechanisms are consistent 
with a given mass action law. We shall show how the rule construction determines the level of coarse 
graining at which the system is viewed. 
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VI. BISTABLE SYSTEMS AND FLUCTUATIONS 



A. The Schlogl model 

Bistable chemical systems show a number of characteristic types of wave propagation processes and it 
is interesting to examine this relatively simple class of systems from the perspective of reactive lattice- 
gas models, not only to illustrate the methods used to construct the automaton rules but also because 
fluctuations can have important effects on the dynamics. 

One of the best-known examples of a model reaction scheme that gives rise to bistable states is the 
Schlogl model. [17] The chemical mechanism of this model is 

A ^ X , 

k-i 

2X + B 1^ 3X . (173) 

k-2 

Letting the concentrations of A, X and 5 be pi, p2 and ps, respectively, the corresponding mass action 
rate law is 

dpi , , , 

— = -kipi + fc_i/92 , 

at 

— = kipi - fc_i/92 + k2P2P3 - K-2P2 , 

^ = -k2plp3 + k.2pl ■ (174) 

We shall use the general formalism of Sec. V to construct several variants of P for the Schlogl model that 
correspond to different coarse grainings of the physical system under different constraint conditions. 



1. Closed system 



Consider the case where all reactive chemical species are treated on an equal footing and only the 
solvent molecules are taken to be uniformly distributed over the lattice. We suppose the system is closed 
so that there are no flows of reagents into or out of the system. The system will relax to a spatially 
homogeneous equilibrium state and no complex behavior is observed. Nevertheless, this case serves as a 
simple illustration of the methods used to construct the automaton rule. 

Given the mass action rate law (174) and using (169) it follows that 

Pi{ol) = -rf{a) + rj"(a) 

P2{a) = r+(a) - rf(a) + r+(a) - r^(a) 

P3(«) = -r2+(«) + r2-(«) . (175) 

where 

rf(<y.) = hkiai , rtM = hk2- 7Ttt2(«2 - 1)«3 , 

(m — ij 

*> 

in 

r^(a) = hk.ia2 , (a) = /iA;_2^— — — ^a2(a2 - l)(a2 - 2) . (176) 

Note that Pt{ol) = since the total number of particles is conserved in the mass action law. At this 
stage, additional considerations are needed to completely specify the reaction probability matrix P. 

If the automaton node is taken to correspond to a very small fluid volume element and the automaton 
time unit corresponds to a small real time interval, the automaton dynamics can be chosen to mimic the 
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individual reactive collision events in the chemical mechanism. Then, the only non-zero elements that 
appear in P are those that correspond to reactions in the mechanism (173). This restriction, along with 
the mean field constraints, allows one to completely specify P. We have for a ^ (i 

P{a,f3) = 6i3-^ai + ll^l32a2-lh3a3iri (ot) - r j*" (a)(5„2m ) ( 1 - l^aim) 

+ f^l3iai-ll^l32a2 + lh3a3irtioi) " J'f (")'^ai™ ) ( 1 ~ l^ajm) 

+ l5/3iai-ll5/32a2'^/33a3(''l'(") " J'f ("))'^ai™ '^"2™ 

+ '^/3i«i'^/32«2 + l'^/33«3-l(''2"(") - r2 (a)6a^rn)(l " 6a2m) 

+ hiaih2a2-lh3a3 + l(r2 (o^) " ''^ («)'5a2™ ) ( 1 " ^a3m) 

+ ^I3iai^l32a2^l33a3-i 



where is the Kronecker delta and 

P(«,«) = l-^^^^P(«,/3) . (178) 

For the input state a = (m, m, m) the condition Pj{{m,m,m)) < for r = 1,2,3 leads to ki = k-i 
and k2 = k-2- Furthermore, for the input states, listed in the argument, the following conditions 
Pi((m, m, as)) = P2((m, m, m)) = Paiiai, m, m)) = imply for the output states /3 that 

_P((m, m, a3),/3) = 0, Pi ^ m , 
_P((m, m, m), /3) = 0, /32 7^ m , 

P{{ai,m,m),f3) = 0, ^3 ^ m . (179) 

Since this scheme allows only those reactions that appear in the mechansim (173) (plus those needed to 
obtain the mean field rate law because of exclusion), this level of description is perhaps the closest that 
the automaton dynamics can be made to correspond to the true microscopic reactive dynamics since the 
fictitious automaton particles undergo reactive collisions like those of the real chemical species. 

If one imagines that the automaton elementary cell corresponds to a larger fiuid volume element or 
the automaton time corresponds to a longer real time interval many reactions can occur in the volume 
element during one automaton time step. In this circumstance the details of the mechanism will not 
play as important a role and the overall kinetics embodied in the rate law will suffice to determine the 
local reactive events in the automaton. It is now a simple matter to write an expression for P that is 
consistent with the overall kinetics and the constraint conditions. We have for a 7^ /3 

3 

P(a,/3) = J2 [p^(«)V«.+i («)V«.-i] n ' (180) 

r = l KjiT 

where 



p+(a) = g+(a)(l - Sa^rn) , 



(181) 



qf{a) = rf{a) , 

g±(a) = rf{a) + r^{a) , 

qt{a) = r^{a). (182) 
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2. Open system 



The reaction dynamics of open chemical systems is far richer than that of closed systems close to 
equilibrium. It is not difficult to extend the automaton dynamics to this situation. Suppose the reaction 
is forced out of equilibrium by flows of one or more of the reagents into the system through the boundaries. 
In the Schlogl model we suppose A and B are fed into the system from reservoirs of these reagents. In 
this case such flow terms must be appended to (174): 

^ = kAiPi - pi) - kipi + k.ip'2 , 

dp2 , , , , 2 ,3 

— = kipi - fc_i/92 + k2P2P3 - k-2P2 , 

^ = -k2plp3 + k.2pl + kB(A-P3) , (183) 

where a superscript is used to denote a constant reservior concentration. This is the usual form of the 
rate law for a reaction in a continuously-stirred tank reactor (CSTR). These flow terms can be represented 
in the mechanism by adding reactions of the form^'* 

kA 

Ao^A, 

kA 

5o ^ 5 , (184) 

ks 

to the Schlogl mechanism (173). If there is one reservoir containg the A and B species then the feed rates 
are equal and = ks = k f . If we let Pj {a) denote Pt{ol) with feed terms then 

P/(«) = P2(«) , 

Pii(x) = P3i(x) + kBip'i-a3) ■ (185) 

The reaction probability matrix may now be constructed for this modifled reaction mechanism and one 
of its possible forms is 

PHa,l3) = P(a,l3) 

+ 6i3^ai + l6l3^a2h3a3hkA 

+ 6j3^ai-l6j32a2^l33a3f^kA{ai — Poi^ami) 
+ ^Piai^P2a2^P3a3 + l^kBP03{^ — l^aam) 

hkB{o:3 - poaSa^m) ■ (186) 

Note that in this model of the far-from-equilibrium system all three chemical species, A, X and B 
fluctuate. 



3. One-variable model 

An even simpler description of the reaction is possible by assuming that the constrained A and B 
species are uniformly distributed over the lattice and their numbers do not fluctuate. This leads to a 



We assume the reactions involving A and B that produce species X are catalysed and do not occur in the 
reservoirs. 
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one-variable description of the reaction dynamics. The results for this case may be obtained from those 
given above in the following way: on the r.h.s. of (176) ai is replaced by pi and as is replaced by 
where pi and are now constant concentrations in the system. In (175) only P'jia'j) is relevant with a 
replaced by a2. Given this value of P2(a2) the reaction probability matrix P(a2, (32) is easily constructed. 
For example, a form for P that includes only those reactions that appear in the mechanism (and taking 
into account exclusion) is for a2 ^ (32 

P{a2,(32) = [(J'r("2) + r^(a2))i5a2-i/32 + {rt{a2) + rj(a2))l5a2+l/3j(l - ba^^m) 

+ [-rt{ci2) + r^{ci2) + r^{ci2) - r^{a2)\^a2m , (187) 

which is the analog of (177). Naturally, it is possible to construct other versions of P{a2,(32) for this 
one-variable case. Now, of course, only fluctuations in the chemical intermediate X are possible. 



4. Simulation results 



One-vanable model 

We begin our discussion of the LGA simulation results by considering the simple one- variable model for 
the Schlogl reaction and use (187) for the reaction probability matrix. [32] In earlier studies a different 
reaction probability matrix was used. [16] We dispense with the subscript 2 for the species X since this 
is the only species whose dynamics is followed. The reaction-diffusion equation is 

dpir,t) _ J J , J , 2/„ ,A „3/„ ,\ , nr72 



kia - k_ip{r, t) + k2bp'{r, t) - k_2P (r, t) + £)V>(r, t) 

(188) 



dt 



6pir,t) ' 

where a = pi and 6 = pg. In the second equality above we introduced the local free energy functional 



^=-J drl^-Vip) + r^D\Vp\ 



(189) 



where the potential V(p) is defined as 



V(p) = -k,ap + '^p' - fp^ + ^p^ . (190) 

Equation (188) is just the time-dependent Ginzburg-Landau equation for a non-conserved order parameter 
field which has been studied often in other contexts. [33] If a random noise term is appended to (188) we 
obtain Model A of critical phenomena [34]: 

%^ = -^^^f^ + €(M), (191) 
at cp(r,t) 

where ^(r,t) is a Gaussian white noise process, 

{ar,mr',t')) = 26{r-r')6{t-t') . (192) 

These considerations place the Schlogl model investigation within a class of well-studied models for both 
the deterministic and stochastic dynamics. Some of the main features of interest in the investigation of 
these systems are: («) the form of the interface separating the stable states, which is known analytically 
for the deterministic model in view of the fact that the dynamics derives from a one-dimensional potential; 
(««) the scaling of the dynamic structure factor on long distance and time scales and its space and time 
dependence on short distance and time scales; and (m) nucleation and growth dynamics. Some of these 
features depend crucially on the existence of fiuctuations; for example, nucleation-induced growth and 
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the short distance and time dynamics, while others, such as the behavior in the scaling regime, do not. 
Thus, the model allows one to determine the ability of the automaton dynamics to accurately reproduce 
the behavior in the deterministic macroscopic regime and to explore the effects of fluctuations on certain 
aspects of the dynamics. 

The automaton calculations were carried out on square N x N lattices with N = 512 and periodic 
boundary conditions using (187) for the reaction probability matrix. The two checkerboard sublattices 
were treated separately in the data analysis. The length scale of the local inhomogeneities in the automa- 
ton simulations is determined by the relative magnitudes of diffusion and reaction. In this study we have 
taken (cf. (160)) £ = 6 (D = 3/2 in automaton units of lattice units squared per time step). This choice 
eliminates most small-scale, short-time reactive recollision events which can lead to significant correlation 
corrections to the steady state average concentrations. [16] 

The homogeneous steady states of (188) are given by the solutions of kia — k-ip -\- k'jhp^ — k-2P^ = 0, 
and are sketched in Fig. 4 as a function of k-i for fixed values of the other parameters (cf. Fig. 4 caption). 



FIG. 4. Schlogl steady states as a function of k-i . The solid line is the deterministic solution while the points are 
the automaton simulation results. The system parameters are: kia = 0.001, 8k2b/3 = 0.095 and 16A;_2 = 0.245. 

The two stable steady states will be called pi and p2 while the unstable steady state will be called Pq. 
Since the dynamics derives from the potential V(p) the stability of the coexisting states can be deduced 
from the relative depths of the potential minima. While the less-stable state is metastable, its lifetime can 
be very long and fiuctuation-induced transitions will occur rarely in simulations on finite systems. The 
deterministic steady state solutions are compared with finite-duration automaton simulations in Fig. 4. 
The initial condition was taken to be a random distribution of particles over the nodes of the lattice with 
mean concentrations corresponding to the deterministic fixed point values. The mean concentration was 
determined from an average over the lattice and over a time T = 2000 steps, following a transient period 
of 1000 time steps. The figure shows that the average concentration corresponds closely to that of the 
deterministic system, with small deviations in the interiors of the bistable domain and larger deviations 
near its boundaries which are due to local fiuctuations in the X concentration. On the time scale of the 
simulations transitions between steady state branches occur only at points very close to the edges of the 
bistable domain. 

Chemical consumption fronts can be generated by taking initial conditions where each half of the lattice 
has average concentration equal to one of the deterministic stable states. The wave front connecting these 
two stable states will move in such a way that the more-stable phase consumes the less-stable phase. The 
point of zero wave velocity corresponds to equistability of the two phases. [17,35] The automaton wave 
front connecting two steady states with the same stability is shown in Fig. 5. 
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FIG. 5. Interfacial profile separating two equally-stable states for k-i = 0.0195. The wave velocity is zero. 
The profile was obtained by averaging over the j/-coordinate on the lattice and is plotted as a function of the x 
coordinate. The solid line is the deterministic prediction. 

The interfacial profile follows easily from the solution of (188) and for the above equistability case it is 
given by 

p(x) = On -\ ^= tanh ^ — , (193) 

where pi = k2b/3k_2, and Q = {{k2bf/3k_2 - 'fe-i}^''^- From Fig. 5 one sees that the automaton 
simulations agree well with the predictions of the deterministic model. In fact, it is also possible to 
derive an expression for the interfacial profile when space is considered as a discrete variable, a situation 
that is more directly applicable to the automaton simulations. In this case the interfacial profile is 
determined by a conservative two-dimensional mapping and points on the interface are given by alternate 
intersections of the stable and unstable manifolds connecting the hyperbolic fixed points that correspond 
to the temporally stable steady states. [36] Concerning the discrete nature of the automaton, one should 
also keep in mind that the behavior of a very sharp front can be anisotropic because the isotropy may 
not yet been recovered at the lenght scale of the width of the front. 

It is also interesting to study the analog of a critical quench where the system evolves from the unstable 
state. The initial decay from the unstable state is strongly infiuenced by fiuctuations; once well-defined 
domains of the two stable phases have formed a deterministic description should be appropriate. In this 
long-time regime the dynamics is driven by the curvature of the boundaries separating the stable phases. 
In an infinite system, although the average order parameter (cf) = p — pg) is zero, domains of arbitrarily 
large size exist in the system. In finite systems different realizations of the evolution process lead to pure 
pI or p2 phases (or mixtures of these phases separated by planar interfaces) but averaged over realizations 
the order parameter will again be zero. The evolution of the system during phase separation may be 
characterized by the nonequilibrium correlation function 

T 

C7(r, t) = {N-'T-' J2 E '^(^' + ^' ^' + ^)'^(^'' ^')) , (194) 

t' = l r'eC 

whose space Fourier transform is the intermediate scattering function S(k,t). Here the angle bracket 
signifies an average over different realizations of the evolution process, and (t)(r,t) is now X]i=i'?2(^iO~Po- 
If the domain size R(t) is the only characteristic length in the system and its evolution is governed by 
interfacial curvature, R{t) ~ t^^^, the intermediate scattering function will satisfy the scaling equation 
[33] 

S{k,t) r~.tF{kt^^^) , (195) 

in two dimensions with F a universal scaling function. 

An example of phase separation following a critical quench is shown in Fig. 6. In the simulation 
k-i = 0.0195 (the equistable point) and the system was uniformly seeded with X particles with average 
concentration corresponding to the deterministic unstable steady state. Sharp boundaries form as time 
increases and slowly deform due to diffusive motion of the interfaces. The final panel in Fig. 6 shows stable 
phases separated by a nearly planar interface. The intermidiate scattering function has been computed 
using reactive lattice gas dynamics and the scaling relation (195) was found to be verified for the late stage 
dynamics. [32] Since this reactive dynamics has a non-conserved order parameter one expects R{t) ~ t^^^ 
scaling. This should be contrasted with conserved order parameter growth where the interface dynamics 
couples to that of the bulk and gives rise to a t^^^ growth law. [33] 
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FIG. 6. One realization of the evolution from the unstable state following a critical quench. 



Three-vanable model 

We next briefly consider the effects of treating the reactive dynamics of all three chemical species in the 
Schlogl model. [37] Experimental studies of pattern formation processes in chemical systems are usually 
carried out under well-controlled conditions using continuously-fed-unstirred reactors (CFUR). [38] In 
these reactors material is fed into the reactor (usually containing a gel or other medium which is designed 
to suppress fluid flow) by well-stirred reservoirs of chemicals whose concentrations are constant. One 
may then envisage a quasi two-dimensional reacting medium, typically a thin layer of a gel or porous 
medium, whose surfaces are in contact with reservoirs containing uniformly distributed chemicals with 
constant concentrations. The reaction probability matrices in Sec. VI A 2 were constructed to treat such 
situations. 

The steady state bifurcation structure depends on the flow rates kA and ks in (184). The steady states 
of (183) are given by the solutions of the equations, 

* fc-iP2 +kAPi .,„.^ 



P3 

where p2 can be obtained from the roots of 



k_2pf + kBpl 



k2P2 + kl 



- (k_ik2kA + kik_2kB + k_2kAkB)pf + (197) 
{kik2kAPi + k2kB(kA + ki)pl) p^^ - (k^ikAkB) p*2 + kikAkBPi = 0. 

In the automaton model ki = k-i and k2 = k-2- For kA = = kf, the steady state bifurcation 
structure is shown in Fig. 7 as a function of ki for two values of the feed rate constant kf. 



FIG. 7. Three-variable Schlogl steady states (c^ = pr/m) as a function of ki for kf = 0.2 and kf = 0.02. The 
solid line is the deterministic solution while the points are the automaton simulation results. The light solid line 
represents the one-variable model results. The remaining system parameters are: 8k2 = 0.0153, /9? = 0.666 and 
p'i = 2.32. 
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For low feed rates one sees that the bistability regime is significantly altered; it is fiattened, shifted 
to higher ki values and spans a larger range of ki. The three- variable automaton results, indicated by 
the points in these figures, are in close accord with the deterministic steady state structure for all feed 
rates studied. The simulation results were obtained by a space average over a 128 x 128 triangular lattice 
and a time average after a transient period. For very high feed rates, kf = 1.0, the one-variable and 
three-variable results are indistinguishable. The hysteresis loops for the A and B species are compressed 
to a nearly line-like form and closely approximate the constant values in the one-variable model. 

Further insight into the character of the dynamics can be obtained by examining the interface that sep- 
arates two stable states in the spatially-distributed system. A planar interface was constructed as above 
by seeding each half of the lattice with particles (now A, B and X particles) whose average concentrations 
corresponded to the stable steady states. Two examples of interfacial profiles at equistability are shown 
in Fig. 8, one for an intermediate feed rate {kj = 0.2) and the other for a low feed rate {kj = 0.02). 



FIG. 8. Three- variable Schlogl interfacial profiles at the equistable point for two feed rates, (a) kf = 0.2 and 
(b) kf = 0.02. The solid line is computed from the deterministic one- variable model and the fluctuation lines are 
the automaton simulation results. 

In the first case the interfacial profile for the X species closely approximates that from the deterministic 
one-variable model (193) and the A and B concentrations do not show noticeable systematic variations 
across the profile. However, for low feed rates no comparison with the one- variable model is possible since 
in the parameter range where the three- variable model exhibits equistability there are no real stable roots 
for the one- variable model. One can see that all three species display well-developed interfacial structure, 
indicating a complete breakdown of the one- variable model. The domain growth and scaling structure for 
early times also show differences when compared with the one- variable model that suppresses fiuctuations 
in the A and B species. [37] These results show that the automaton dynamics is able to accurately describe 
the fiuctuating dynamics of bistable chemical systems and, in regimes where fiuctuations are unimportant, 
automaton results are in accord with deterministic predictions. 
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B. Fluctuations and second moment constraints 



In previous sections we pointed out that lattice gas automata provide a microscopic approach to cooper- 
ative phenomena and that by construction the automaton possesses intrinsically spontaneous fluctuations. 
Therefore the LGA method should be well suited to investigate the role of fluctuations on the macroscopic 
behavior of complex systems. This indeed is an important point at least for two reasons: (i) reaction- 
diffusion systems exhibit macroscopic and mesoscopic space- and time-dependent behavior which often 
is triggered or influenced by fluctuations; (ii) these fluctuations are usually not easily accessible to other 
methods. Now it is crucial that the simpliflcation introduced in the microscopic LGA description should 
manifest neither at the macroscopic level nor at smaller (mesoscopic) scales where the fluctuations can 
play an important role. So in order to obtain relevant information, we should have some guarantee that 
LGA fluctuations capture the essential aspects of actual fluctuations . This problem can be addressed 
at least within the limits of the linear theory of reaction-diffusion systems [39]. For simplicity we consider 
a one- variable system. 

A standard approach for non-equilibrium steady states (far from a bifurcation point) is the Landau 
method where a process-independent noise is added to the phenomenonlogical RD equation linearized 
around a steady state {p =p* -\- <j)) 



dt 



-nc^iv, t) + OV'^c^iv, t) + ez)(r, t) + ^r{v, t), 



(198) 



with K the linearized reaction rate coefficient and where the diffusive noise and the reactive noise 
are taken to be additive and whose mean values and correlations are given by 



< CD(r,t) > 

< CR(r,t) > 

< CD(r,t)CR(r,t) > 
<^Dir,t)^Dir',t') > 
<Uir,tKRir',t') > 









-AD6{t-t')V^6{r-r'), 
AR6(t-t')6(r-r'). 



(199) 



Now the main problem is to characterize the statistical properties of (f). In general, these properties 
depend on the initial condition which may be stochastic; however, in the long time regime all details 
of the initial condition fade away and it becomes possible to characterize (f) without reference to any 
particular initial state. In this regime - the only regime considered here - the flrst moment of (t)(r,t) is 
zero: 



< <j>(r,t) >= 



(200) 



which follows from (199) and the linearity of (198). A less trivial characterization of (f) is given by the 
density fluctuation correlation function 



Gir,t)={4>ir',t')4>ir' + r,t' + t)) 



(201) 



Here we are mainly interested in the static correlation function G(r) = G(r, t = 0) which can be evaluated 
from (198) and (199): 



^^For instance the thermal lattice gas model constructed by Grosfils, Boon and Lallemand [40] has been shown 
to exhibit spontaneous fluctuations whose correlation function is in agreement with the dynamic structure factor 
measured in real fluids. 
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where the function Hd(z) depends on the space dimension d. For d = 1, 2, 3 we have 

H2{z)= —Kq{z) (modified Bessel function); 

Mz)= \i^f""^ (203) 

respectively [39]. In the automaton, the correlation function C{v,t) is defined as in (201): 

C{v,t) = {^{v',t')^{v' + v,t' + t)) , (204) 

but here r and t take discrete values and <f)(r, t) = '}2'i=iVi{^ > Po- In ^ simulation, the static correlation 
function C(r) = C(r, 0) can be evaluated by space and time averaging 

T+t 

C{r) = N-'T-' J2 E '^(^' + ^' ^') , (205) 

t'=t r'eC 

where t must be large to reach the long time regime, and where T must also be large to have a good 
estimate of C; averaging over different realizations can also be performed. 

In order to compare (205) with the predictions of the Landau approach, we need to evaluate the 
factors Ao and Ar appearing in (199). To obtain the value of Ao , we consider a node in an equilibrium 
configuration at the steady state density p* , and we evaluate the variance of the number of particles on 
that node 

m 

Var[^ryi] = p* {1 - p* /m) . (206) 

8 = 1 

Then, we evaluate a similar variance in the Landau approach by integrating G(r) over a domain of volume 
V corresponding to a node of the lattice. In the limit of infrequent reactions k ^ 0, we obtain 

Yai [p] = V Ad /2D . (207) 

Indentifying the two variances (206) and (207), we have 

AD=2yp*(l-p*/m) . (208) 

The value of the factor Ar depends on the probability matrix P: a small value is obtained when the 
changes in particle number (/3 — a) are concentrated around their average value; conversely Ar has a large 
value when the changes (/3 — a) are dispersed around their average value. Now the phenomenological rate 

F(P) = J2(P-a) (p/^r (1 - p/mr-"P(a(3) , (209) 

is the average change in particle number produced by the action of the operator C; analogously Ar is 
given by the variance of the change in particle number ^® 



The variance is calculated with respect to the grand canonical ensemble distribution. The microcanonical 
distribution gives the same result in the strong-diffusion limit [22]. 
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Ar 



a ,13 




(210) 



This is an important relation which allows one to control the level of the reactive noise in the automaton. 
Indeed, for a given macroscopic rate F(p), it is possible to independently prescribe the function Ar(p), 
at least within some limits [39]. 

Figure 9 shows the static density correlation function G(r) for different microscopic dynamics cor- 
responding to a fixed linear macroscopic rate law. Notice that the system exhibits negative spatial 
correlations when the reactive noise amplitude is minimized (Ark~^ < A£)D~^) and positive correla- 
tions in the opposite case. The agreement between simulation data and theoretical predictions is seen to 
be excellent. The effect of dimensionality is also shown in Fig. 9: clearly the range of spatial correlations 
decays with the dimension of the system as intuitively expected. 



FIG. 9. Left panel: Spatial fluctuations correlation function G{r) in a two-dimensional reactive lattice gas; 
parameters : Ad / D = 0.1 and Ar/ k = 0.05,0.1, and 0.225 (see text for explanation). Error bars are smaller 
than the size of the symbols. Right panel: Static density fluctuations correlation function G(r)/G(l) in one- 
(triangles), two- (circles) and three- (crosses) dimensional lattice gases. Symbols are simulation data and solid 
curves are obtained from (202) and (203) . 

These observations show that the steady state fiuctuations in reactive LGA are consistent with the 
predictions obtained from the Landau theory (at least within the limits of linear theory). This is an 
important result in that two very different approaches yield the same results. 
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VII. EXCITABLE MEDIA AND SPIRAL WAVES 



Chemical waves in excitable media are some of the most important types of wave processes seen in 
physical systems because they play a major role in the operation of the nervous system, the heart as well 
as a variety of other biological processes. [42] In addition, they have been shown to be involved in cat- 
alytic oxidation processes on metal surfaces [43] and other chemical systems. The Belousov-Zhabotinsky 
reaction [44] can be carried out under excitable conditions and typically displays such chemical waves. 
[45,46] 

A chemical system is excitable if it has a stable fixed point and responds to perturbations in the following 
way: initial conditions obtained from small perturbations of the fixed point give rise to trajectories that 
make small excursions in phase space or return directly to the resting state in a short time. Perturbations 
that exceed a threshold value give rise to trajectories that make a large excursion in phase space before 
returning to the resting state. During such long excursions the system is refractory and insensitive to 
perturbation. If the system has spatial extent, the excitable medium will respond to local perturbations 
by producing waves of excitation with various geometries that travel through the system. The spatial 
characteristics of the wave refiect the temporal dynamics of excitability described above: the wave front 
corresponds to the system in its excited state while the tail of the wave is refractory as it recovers to 
the fixed point state. As noted above many chemical systems show such excitable behavior and we shall 
focus on one such system below: the Selkov model [47]. 

Our aim in this subsection is to give a simple illustration of the fact that reactive lattice gas automata 
can be used to study excitable media and to examine how concentration fiuctuations can give rise to the 
spontaneous formation of waves in excitable media and infiuence their propagation. 



A. The Selkov Model 

The Selkov reaction [47] 

A ^ X , 

k-i 

X + 2Y ^ 3Y, 

k-2 

Y^B, (211) 

k-3 

was originally constructed to model certain aspects of glycolysis. Here we are not concerned with its 
basis in real chemistry but simply use it to illustrate wave propagation in excitable media. The mass 
action rate law of the above reaction is 

kia - k-ipi - k'2Pip\ + k-'2p\ , 

k-sh - k3p2 + k2Pipl - k_2pt , (212) 

where pi and p2 refer to the concentrations of the intermediates X and Y, respectively. The reversible 
version of the reaction was investigated in the well-stirred limit by Rehmus et al. [48] and was shown to 
display a wide variety of different attractors, ranging from fixed points of various types to bistable states 
and limit cycles. We assume that the concentrations pj^ = a and pg = h oi the A and B species are held 
fixed by external constraints and are taken to be the control parameters along with the rate constants 
in the model. They are treated in the same mean field approximation as the A and B species in the 
one-variable Schlogl model. It is possible to choose these parameters such that the system is excitable. 
The excitability of the dynamics is demonstrated in Fig. 10 which shows the evolution from two initial 
states that are above the threshold. 



dpi 
dt 

dp2 
dt 
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FIG. 10. Two phase plane trajectories for the Selkov model in the excitable region. Both trajectories start 
from initial conditions (open circles) that exceed the threshold for excitability and undergo long excursions in 
phase space before return to the fixed point, which is indicated by a heavy dot. The concentrations are scaled by 
(^2/^3)"'^^ and the time is scaled by ^3. 

Note that the system does indeed make long excursions in phase space before return to the fixed point. 
Initial conditions that lie close to the fixed point relax directly back to it and do not make such large 
phase space excursions. 

One can apply the general formalism of Sec. V to construct a variety of reaction probability matrices 
for the reactive dynamics, as illustrated above for the Schlogl model. From the mass action rate law 
(212) and (169) we obtain the relations 

Pi{a) = rf{a) - rj"(a) - rj(a) + r^(a) , 

P2(«) = r+{<y) - r2-(«) + r+{<y) - (a) , (213) 

where 

(oi) = hkia , ^1 i*^) = hk-iai , 

r+(«) = ^aia2(a2 - 1) , i<y) = ^2(^2 - l)(a2 - 2) , (214) 

r||"(a) = hk-^h , ''3 (<^) = hk2,a2 ■ 

If we construct P{oL,(j) so that only the steps in the Selkov mechanism are allowed (subject to exclusion 
constraints) we find for a 7^ /3 [18] 

_P(a,/3) = (^/J^a^l [r^{a)6p^ai + l + J'f (")'^/3iai-l] (1 ~ ^aim) 

+ (rf(a) - r'l{oL))8a^r 

+ I [''2"(")'^/3i«i-l'^/32«2 + l + i0L)hiCi + lh2<y2-l\ (1 - '5aim)(l " ^a^m) 

+ (rj(a) - r^(a))((5„im - 6a^rn) 

+ I [''3 (")'^/32«2-l + f"-!" (")'^/32«2 + l] (1 - ^a2m) 

+ (''3 (") - rt{a))6a^rn^ 

(r2+(«)-r2-(«)) . (215) 

If instead we simply demand that the mass action rate law be satisfied but the individual steps in the 
mechanism need not be taken into account, one of the possible choices for the non-diagonal elements 
(a ^ f3) oi P{a,f3) is 

2 

] n ' (216) 

T = l k(k^t) 



where 



pt(a) = qt(a)(l - 8a^m) , 

= - ) (217) 
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qf(a) 



rf{a) + rj(a) , 
rtio") + ''3 (") • 



(218) 



B 



Simulation results 



Simulations of the Selkov excitable medium were carried out using (215) on square lattices with m = 4 
with the 4-equi-rotation rule to scramble the velocities. [18] Kinetic parameters were chosen to lie in the 
excitable region and all diffusion coefficients were selected to be equal. Their magnitudes were controlled 
by varying £1 and £2 which were taken to be equal (see Automaton rule in section V). In the simulations 
described below the system was first allowed to relax to the fixed point. If diffusion is sufficiently strong 
the system will be spatially homogeneous on average, exhibiting only small local fiuctuations away from 
the fixed point value. A growing ring of excitation can be initiated in the automaton by selecting the 
average concentration within a disk of radius R to differ from the steady state concentration. Provided 
the concentration differs sufficiently from that of the steady state and the radius is larger than some 
critical value Rc, this local concentration perturbation will excite nodes in the perimeter of the disk 
leading to a large concentration change before relaxation back to the fixed point concentration occurs. 
This process will generate a traveling wave of excitation that moves out from the disk. The chemical 
wave has a refractive tail and leaves behind unexcited medium. An example of this ring growth is shown 
in Fig. 11. 



FIG. 11. Lattice gas automaton growth of a ring of excitation for the Selkov modeL The system parameters 
are: ha = 0.00008385, k-i = O.lfcs, ^3 = 0.0005, k2 = k-2 = 0.01 and k-ab = 0.000002326 with Di = D2 = 3/4. 
Note that the values of the fcr's are indicative of the sensitivity of the system behavior to the parameter values. 
The simulations were carried on a 1024 x 1024 square lattice. 

Notice that apart from small fiuctuations, the rings are circular in shape, while simple cellular automa- 
ton rules often produce square waves with sharp corners. [49,50] 

At this point is worth stressing that the mass action rate law follows by construction from the automaton 
rules. As a result all wave characteristics like wave dispersion and velocity are automatically determined 
once the mass action rate law, or better, the reaction mechanism is specified. This makes the construction 
of lattice gas models that refiect the true reactive dynamics and wave propagation processes a simple and 
straightforward task for any reactive system. Complex automaton rules that are designed to incorporate 
such features have been constructed to mimic the details of excitable kinetics. [51,52] 

In the above simulation the spontaneous fiuctuations have a spatial correlation length that is small 
compared to the critical radius R^. However, if we reduce D, the critical radius - which scales as R^ ~ 
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D/vc ~ D^l"^ - will decrease until spontaneous fluctuations can produce perturbations that exceed R^. 
Here is the velocity of a planar chemical wave, which depends on the diffusion coefficient as ~ D^l"^ . 
[53] Once fluctuations can produce threshold-exceeding perturbations, the medium will spontaneously 
generate waves of excitation. This is shown in Fig. 12. The various frames in this flgure show how 
fluctuations can continuously produce waves of excitation that spread over the surface of the 2-d system. 



FIG. 12. Spontaneous nucleation of rings of excitation. Parameters are the same as Fig. 11 except that 
Di=D2 = 1/4. 

The generic form for a chemical wave in an excitable medium is a spiral since a planar or circular 
wave front fragmenting its free ends will curl to form the core of a spiral wave, provided the medium 
is sufficiently excitable. There is an extensive literature on the conditions for spiral wave formation in 
excitable media. [54] In addition, the dynamics of the spiral core itself need not be simple and can undergo 
bifurcations to periodic, quasiperiodic and perhaps chaotic patterns. [46] In three-dimensional media the 
spiral core, a topological defect, is drawn out into a filament and 3-d spirals such as scroll waves are 
formed. [54] Reactive lattice gas automata allow one to investigate such phenomena as well how they are 
affected by fluctuations. While detailed, quantitative investigations of the effects of fluctuations on spiral 
waves in excitable media have not been carried out, preliminary studies indicate some of the phenomena 
to be expected, such as for example, the irregular wave fronts seen in the simulations described below. 

Here we simply demonstrate that the reactive lattice gas automaton can yield spiral waves with the 
correct qualitative features. A spiral wave can be generated from a ring of excitation by shearing the ring 
to produce free ends which can form the cores of two counter rotating spiral waves. 

Figure 13 shows the formation and evolution of a spiral wave pair. 
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FIG. 13. Formation and evolution of a pair of counter rotating spiral waves produced by removing one half of 
a growing ring. Parameters are the same as in Fig. 11. 

In the simulation the chemical wave was sheared by simply randomly reseeding one half of the lattice 
with particles whose average concentration corresponds to the steady state. This produces a semi-circular 
ring with two free ends. As expected the spiral wave pair continually regenerates itself after the wave 
fronts collide. 
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VIII. TURING PATTERNS 



In a Turing bifurcation [55] a homogeneous steady state loses its stability and an inhomogeneous state 
is formed. While these bifurcations are believed to be relevant in many biological pattern formation 
processes [53], it is only recently that Turing pattern formation has been observed in gel reactors [56,57]. 
These experimental observations have led to increased research activity on both Turing pattern formation 
and more complex scenarios involving interactions between Turing and other bifurcations. [58,59] It is now 
widely believed that immobilization of some species by the gel is responsible for the diffusion coefficient 
differences that drive the Turing instability. [60,61] One of the characteristic features of Turing patterns 
is that a regular or quasi-regular structure with a given (set of) wavelength(s) can develop in unconfined 
systems with no imposed characteristic macroscopic length (in contrast to most typical situations in 
hydrodynamic instabilities). The wavelengh of the Turing pattern is an intrinsic property of the system: 
it depends on the kinetic parameters and on diffusion coefficients. 

Turing bifurcations are usually described in terms of the activator-inhibitor kinetics of a two- variable 
reaction-diffusion equation. [53] Suppose p =< pi,P2 > is the vector of concentration variables, F(p) 
describes the kinetics and D is a diagonal diffusion coefficient matrix, Djjt = Djbjjt . The conditions 
for a Turing bifurcation are as follows. Suppose A = {dY / dp)p=p* is the Jacobian matrix evaluated 
at the steady state p*, and B = A — fc^D is the matrix that governs the linearized evolution of the 
Fourier transform of the concentration fields. If An > and A22 < then species X = Xi is the 
activator and species Y = X2 is the inhibitor. A Turing bifurcation will occur if detH = 0, TrB > and 
A11D2 + A22D1 > 0. The wavenumber at the bifurcation is 

"^''"^ ' ' (219) 
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Furthermore, since the bifurcation must occur from a stable homogeneous steady state we must have 
D1/D2 < 1, i.e., the diffusion coefficient of the inhibitor is greater than that of the activator. The critical 
diffusion ratio at the bifurcation is 

^ = A^^ (detA - A12A21 + 2(Ai2A2irfetA)i/2) . (220) 
The unstable wavevectors are determined from the dispersion relation <T{k) = ^^(A), where 

\{k) = ^TrB ± ((TrB)2 - Ademf^ . (221) 



A. Selkov model simulations 



It is not difficult to realize the above conditions for the Selkov model. The system parameters may 
be adjusted to lie in the region where the stable fixed point of the spatially-homogeneous system is 
a focus and is close to the Hopf bifurcation point. For a given set of rate coefficients and control A 
and B concentrations the diffusion coefficient ratio can be adjusted to yield a Turing bifurcation in the 
reaction-diffusion equation. A set of rate constant parameters that achieves this is kia = 0.002656673, 
fcs = 10A;_i = 0.00665, k2 = k_2 = 0.015, and A;_3& = 0.000531334. With this set of parameters the 
critical diffusion coefficient ratio for the onset of a Turing bifurcation is D2/D1 = 16.2. 

The lattice gas simulations of this bifurcation [18] were carried out on a pair of coupled square lattices 
using the 4-equi-rotation rule to scramble the velocities and the reaction probability matrix (215), just as 
in the excitable medium simulations. Only the kinetic parameters in (215) were changed to correspond 
to the Turing region. The diffusion coefficient ratio was changed by earring out different numbers of 
propagation and rotation steps for the Xi and X2 species so that D2/D1 = £2/^1- When the Selkov 
automaton dynamics is carried out with Di = D2 and the above set of rate constants no macroscopic 
pattern formation was observed. The system exhibited small local fluctuations about an average concen- 
tration corresponding to the steady state. When the system was started from a uniform random initial 
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state with average concentration different from tliat of the steady state, the average concentration was 
observed to decay in an oscillatory fashion to the steady state. These simulations confirm the stability of 
the steady state and its focal character. When the same type of simulation was carried out with the dif- 
fusion ratio selected to lie above the predicted Turing instability the system was observed to evolve to an 
inhomogeneous steady state. Figure 14 shows the results of a simulation with Di = 1/2 and D2 = 25/2. 



FIG. 14. Formation and evolution of a Turing pattern. The system size is 1024 x 1024 and the parameter 
values are given in the text. Here Di = 1/2 and D2 = 25/2. The figure shows a realization of the evolution 
corresponding to one of the two checkerboard sublattices. 

Spatial inhomogeneities develop quickly from the random initial state. The system undergoes bulk 
oscillations during the evolution to the inhomogeneous steady state. The final state is a hexagonal 
pattern of spots distorted by molecular fiuctuations. The magnitudes of these local fiuctuations are 
controlled by the system size and the magnitudes of the diffusion coefficients, which determine the local 
diffusive length scales. By maintaining the same diffusion coefficient ratio but reducing the magnitudes 
of both Di and D2 one can effectively destroy the Turing pattern by fiuctuations. The characteristic 
wavelength of the pattern can be made sufficiently small so that the number of "molecules" in a volume 
with characteristic size corresponding to this wavelength is small enough to give rise to fiuctuations that 
destroy the pattern. This is illustrated in Fig. 15. While the system still shows strong inhomogeneities 
in the stationary regime, now the pattern is far more irregular. Thus the reactive lattice-gas automaton 
simulations can be used to probe both system size and fiuctuation effects on Turing pattern formation. 
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FIG. 15. Formation and evolution of a Turing pattern. Same as Fig. 14 except that Di = 1/4 and D2 = 25/4. 

Often autocatalytic reactions take place in geometries with small linear dimensions. Biological cells 
are perhaps the most familiar examples where the cell dimensions can be quite small, of order 0.1 //m. 
In such small geometries fluctuations are likely to be quite important and it is interesting to speculate 
as to whether concentration inhomogeneites can develop within the cell as a result of the operation of a 
Turing or Turing-like mechanism. Lattice-gas methods provide a tool for the study of such problems. [62] 



B. Maginu model simulations 

Another model which exhibits Turing structures is the Maginu model [63]. This mathematical model 
is defined by a set of PDE's that belong to the class of RD equations: 

dxi o 2 

— — = Xl - - X2 + -Dl V Xl , 

ot 

-^ = (x^-kx2)/c + D2Vx2, (222) 

with c > and < A; < 1. It should be noticed that the corresponding rate law 

dxi 3 

— = Xl - Xi/3 - X2 , 

^ = (xi- kx2)/c . (223) 

is not a mass action rate law because each rate dxi/dt (i = 1,2) can be negative when Xj =0. As a 
consequence, the set of equations (222) cannot be simulated directly by a lattice gas automaton. However, 
it is possible to perform a linear transformation on the variables (xi, X2) so that a lattice gas automaton 
can be constructed for the transformed equations. This can be done in various ways, e.g. 

pi = 1/2 + Xl/ ^12(1 + l/k), 

P2 = l/2 + fca;2/yi2(TTT7fc) , (224) 
which yields to a transformed rate law 

^ = -iaipl + Gaipl - 02/91 -I- 03(1 - P2) , 

^ = a.ipi - p2) , (225) 

with 

ai = {l + l/k)/c , 02 = (2 i/k)/c , 03 = l/fcc , 04 = k/c . (226) 
Using the general formalism of Sec. V to construct a reaction probability matrix P we obtain from (225) 

Pi(a)/h = — 4aiai(ai — l)(ai — 2) 



(m — l)(m — 2) 
in 

6aiai(ai - 1)- — - a2ai + 03(1 - 02) , 

(m — 1) 

P2(a)/h = ai(ai - 02) , (227) 

where the scaling factor h has been incorporated. These relations, together with the mean field condi- 
tion (169), do not specify a unique matrix P. Additional requirements about higher moments must be 
provided in order to fully determine the matrix. The following rules are used: 
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1. for any initial configuration a, the changes in particle numbers are performed independentely for 
each species: 



(228) 



2. in any reaction, the number of particles of each species changes at most by one unit: 



Piiai,(]i) = 0, 



if \ai 



A-l > 1 ; 



(229) 



3. for any initial configuration a, creation and anihilation transitions are exclusive within each species: 



The purpose of these conditions is to reduce the number of different reactions occuring on the lattice, and 
to avoid large fiuctuations around the mean field rates. The above rules specify completely the matrix P 
except for the value of the scaling factor h which remains free at this stage. 

Now we turn to the definition of the diffusive part of the dynamics. The goal is to construct a lattice gas 
rule with enough fiexibility to allow for smooth changes in the value of the different diffusion coefficients 
Dr ■ This is important because a small change in the value of a diffusion coefficient can trigger the Turing 
instability. To obtain this fiexibility, we used an automaton rule 



where the strength of the shuffiing performed by each velocity randomization Rr is controlled by a 
corresponding parameter g^. Here, Rj{qj) is the velocity randomization operator (84) with = 1 — 3qr 
and pi = p2 = ps = Qr ■ Using the result (94) and taking into account the value of f^, we find the following 
expression for the diffusion coefficient expressed in automaton units (lattice units squared per time step) 



Large values of ( > 5 in automaton units) are often needed to simulate Turing structures. From (232) 
we see that two parameters (£r and g^) can be used to assign a given value to Dr- A good strategy is 
to choose the largest possible value for £r , and keep the fiexibility provided by the continuous parameter 
Qr for the fine tuning of the diffusion coefficient; in this way the diffusive behavior emerges at a length 
scale of a few lattice units. Indeed, the opposite strategy (£r = 1, and a small value for g^) should be 
avoided because in this case the discrete lattice structure is more likely to appear in the Turing structure 
which will then be different from the prediction of the RD equations. [22] The only way to avoid this 
is to make sure that the typical length scale of the Turing structure is larger than the scale required 
for the emergence of a good level of isotropy in the diffusion. This requires larger lattices and more 
updates. It should be stressed that the parameters {£r : r = 1, . . .n} must have the same parity in order 
to maintain the two checkerboard subsystems uncoupled. Neglecting this requirement usually produces 
strong artifacts. [22] 

In some circumstances, the discreteness of time can also become apparent in a Turing structure obtained 
with an automaton; in particular, the order in which the basic operators Pr, Rt and C are combined to 
construct £ can be important. For instance, if the automaton rule is (231), the system is always observed 
after the diffusive phase of the dynamics and before the reactive transformation; in some sense the Turing 
structure is observed through a diffusive filter. On the other hand, the equaly acceptable automaton rule 



Pi{ai, {cti — 1)) 7^ is not compatible with Pi{ai, {cti + 1)) 7^ . 



(230) 




(231) 




(232) 




(233) 
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produces sharper Turing structures. The differences between the two automaton ruies shouid vanish 
when the densities vary smoothiy over iarge enough space and time scaies. 

The simuiations presented here were carried out on square N x N iattices with N = 256 and periodic 
boundary conditions. The two checkerboard subsystems were aiways considered separateiy (see Sec. III). 
Linear stabiiity anaiysis of the Maginu eq uations shows that the homogeneous steady state can become 
unstable by spatial destabilization when Di/ D'j < (1 — Vl — k)^/c/k for < c < A; [64] in which regime 
the simulations presented here were performed. 

Figure 16 shows a wormlike Turing pattern obtained from the destabilization of a system initially 
prepared in the homogeneous steady state (according to the rate equations (225)). 



FIG. 16. Turing pattern (pi in one of the two checkerboard subsystems) obtained from the destabilization of 
the homogeneous unstable steady state in the Maginu model (225). Lattice size: 256 x 256; k = 0.9, c = 0.45, 
h = 1/40, li = 1, I2 = 9, qi = 1, g2 = 0.6429, (Di = 0.25, D2 = 4.75); state after 14000 time steps. 

Figure 17 shows a perodic Turing pattern obtained from a perodic initial condition. A quantitative 
comparison of this structure with the prediction of the phenomenological RD equations shows very good 
agreement [22]. 



FIG. 17. Turing pattern (pi in one of the two checkerboard subsystems) obtained from a periodic initial 
condition. Lattice size: 256 x 256; k = 0.9, c = 0.45, h = 1/60, li = 1, I2 = 9, qi = 0.4, q2 = 0.21176, (Di = 1, 
D2 = 19). 
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IX. OSCILATIONS AND CHAOS 



A. Limit cycle oscillations 



Under certain conditions, the Maginu rate equations (223) (or equivalently (225)) exhibit oscillatory 
behavior. A set of parameters corresponding to this regime is (k = .9 and c = 2). Simulations of the 
automaton were carried out for these values with Di = D2. The automaton is prepared in a homogeneous 
state and the spatial average of the density of each species is recorded during the system evolution. Phase 
trajectories in the (pi , /92)-plane are obtained in this way (separately for each checkerboard subsystem). 
Figure 18 shows how a trajectory is initially attracted by the limit cycle predicted by the phenomenological 
rate equations, then shrinks progressively to a smaller cycle whose amplitude depends on the value of the 
scaling parameter h while its frequency remains locked to the predicted phenomenological value. When 
h is large, reactive fluctuations are strong, and diffusion can no longer maintain spatial coherence over 
the whole system. As a result, different zones of the lattice become progressively phase shifted with 
respect to each other. As time evolves, the phase decoherence increases up to a point where a balance 
is reached between the source of decoherence (fluctuations) and the homogenization due to diffusion. So 
the observed limit cycle is contracted due to the combined effects (resulting from spatial averaging) of 
many out-of-phase local limit cycles with normal phase space size. 

The values of the frequency of the limit cycle obtained from the power spectrum of the LGA simulation 
data ipiii), p'zii)) are in agreement with the phenomenologically predicted value, Eqs. (222). However 
there is a slight discrepancy for large values of h which can be interpreted as the result of the strong 
gradients and could also be the consequence of insufficient local diffusive equilibrium when the ratio of 
reactive collisions to elastic collisions becomes large. 



FIG. 18. Phase trajectories of limit cycle behavior in the Maginu model. Left panel: LGA simualtions (one of 
the two checkerboard subsystems). The difference between the outer ring and the inner ring shows the shrinking 
effect when the ratio of reactive collisions to elastic collisions is increased, i.e. A = 10 versus h = 2. Lattice size: 
64 X 64 . Right panel: limit cycle obtained by numerical integration of the Maginu model (225). 
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B. Oscillations and chaos 



Chemically reacting systems provide some of the best-characterized examples of deterministic chaos and 
a considerable amount of both experimental and theoretical research has been devoted to an elucidation of 
its properties and the mechanisms responsible for its appearance. [65] In its simplest form deterministic 
chemical chaos is manifested in well-stirred systems with three or more chemical species. All of the 
theoretical investigations of such temporal chaos have been carried out using the ordinary differential 
equations of mass action kinetics. In this section we examine "deterministic" chemical chaos from the 
mesoscopic point of view of the reactive lattice gas automaton. 

This more microscopic picture of the chaotic dynamics will allow us to consider a number of fundamental 
questions concerning these systems. We shall show how it is possible to construct a mesoscopic reactive 
dynamics whose mean field description yields deterministic chaos. [66] On the basis of this dynamics 
we can then investigate how large the system must be before one obtains a recognizable attractor. The 
interplay between system size, diffusion and reaction and how they affect the structure of the attractor 
will be examined. 



1. The Willamowski-Rossler model 



The Willamowski-Rossler model [67] is an example of a chemical mechanism that gives rise to a chaotic 
attractor. Since it is based on a scheme with mass action kinetics it is especially suitable for the investi- 
gation of the microscopic basis of chemical chaos. 

The Willamowski-Rossler reaction mechanism is 

Ai+ X ^ 2X, X + Y ^ 2Y, 

A5 + Y 1^ A2, X + Z 1^ A3, (234) 

k-3 k-4 

ks 

A4 + Z ^ 2Z , 
k-s 

where k = {k±i : i = 1, • • • , 5} is a set of forward and reverse rate constants. In (234) A = {Aj : 
j = 1, • • • , 5} are a set of species whose concentrations are fixed by constraints. The corresponding mass 
action rate law for the concentrations of remaining species X , Y and Z is (t = X,Y, Z = 1,2,3) is 

dpi 2 2 

— = Kipi - K-lPi - K2P1P2 + K-2P2 - K4P1P3 + K_4 , 
dp2 2 

— = K2P1P2 - K-2P2 - K3P2 + K-3 , 

^ = -K4/91/93 -I- K_4 -I- K5/93 - K_5/93 . (235) 

In (235) the constants Ai have been incorporated into the new set of rate coefficients {k} = {K±i : i = 
l,---,5}. 

One can recognize in the first three reactions in (234) a Lotka-Volterra oscillator involving the auto- 
catalytic species X and Y . This Lotka-Volterra oscillator is coupled via X to a "switch" between the 
autocatalytic species Z and X . A suitable coupling between the Lotka-Volterra subsystem and the switch 
element leads the observed oscillations of the full Willamowski-Rossler system. [68] 

The reaction probability matrix corresponding to the Willamowski-Rossler model used in this study is 
[66] 

P(a,f3) = [p^(a)6i3^a^ + i + (a) Sp^a^-i] 8p^a^8p^a^ 

+ [P2"(") '^/32«2 + l +P2 (") '^/32«2-l] '^/3i«i'^/33«3 (236) 
+ [P3 (")'^/33«3 + l +P3 (")'^/33«3-l] hiaih2C,2 ) 
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for OL ^ (j, and 



The explicit expressions for p^{a) are given below: 



P(«,«) = l-^ P(«,/3). (237) 



= q+{oL){l - 6a^rn), Pi («)= Qi («) " qt(a)6a^ 



m 1 



»2"(") = 92"(")(1 - ^a^m), P'j («)= 92 (") " 



= 9j(«)(l - 6a,m), Ps (a)= 93 (") - 9j(«)'5. 



with 



it , mK-2 / i\ , 91 mK_i / i\ , / , \ 

-jL = Kitti + ^ _ ]^ 02(02 - 1) + -j^ = ^ _ ]^ ai(ai - 1) + (K2a2 + K4a3)ai 

99" , 99" niK-'i , 1N , 

= K2aia2 + K-3, = rn - l "2(a2 - 1) + ^3^2 , 



9i , 93 , mK_5 . 

-^=K5a3 + K_4, = K4aia3 + ^ _ ]^ 03(03 - 1) 



When this reaction probability matrix is substituted into the automaton mean field equations one 
obtains the Willamowski-Rossler mass action rate law. Since the full automaton dynamics is not mean 
field, we can now use the reactive automaton to investigate the dynamics of this reacting system from a 
mesoscopic point of view that incorporates fiuctuations. 



2. Spatially homogeneous system 



Before examining the dynamics of the Willamowski-Rossler system using the full automaton dynamics 
it is interesting to first examine a simpler situation where the system remains fully mixed or well stirred 
throughout the evolution. In view of the description given earlier, the particles are binomially distributed 
in the well-stirred limit. In this case a Markov chain equation can be written for the time evolution of 
the probability P(ii,k) that there are n = («i,«2,«3) particles of species Xi on £1, X2 on £2, and Xs 
on £3. The Markov chain reads, 

P(ii,k+l)- P(ii,k)= [W(n\ii')P(ii',k) -W(n'\ii)P(ii,k)] , (238) 

The transition probability 1/F(n'|n) from n to n' can be written as 

3 

P^(n'|n)=n E Prinr+,nr-\pit)) 

r = l T^r+'^r- 

_(nr+-nr-=n'r-"r) 

where _P^(n^_|_, jp) is the probability that there are particle transformations on lattice £^ that 
increase the particle number of species Xr by one and particle transformations on lattice £^ 

that decrease the particle number of species Xr by one. Here we have used the fact that in the 
Willamowski-Rossler mechanism the reactions involve particle number changes of ±1. The explicit form 
for Pr(nr+,nr-\p) is 

Pr(n+,n.\n) = c(n+, n_, W)(P/)"+(P-)"-(l - P+ - p-)^-"+-"- , (240) 
where A/" = N^. The factor 



(239) 
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c(n_|_, n-,M) 



n_)! 



(241) 



takes into account the number of different ways of assigning the transitions ^ + 1 and n_ 
transitions ^ — 1 to possible positions on the lattice. The _P+ and P~ are transition probabilities 
averaged over a binomial distribution on £^ and are given by 



c^(l-c.)^ 



(242) 



One may again obtain the mass action rate by computing the average value of Prik) from the master 
equation (238) in the mean field limit: 



E[Apr] = j^ J2 (n+-n_)Pr(n+,n_\p) = P+ -P- =hRr(p) 



(243) 



n_|.+n_ <Af 



We may also compute the standard deviation of the net particle number change per node on a(Apr), 
in the same approximation and find 



E 



E(Apr) 



Pr(n+,n_\p) 



n_|.+n_ <Af 

l:[p+ + p--(P+-p-f] 



(244) 



This result shows that the fiuctuations scale as Af~^^^ as expected. 

One way to simulate this well-stirred dynamics is to introduce a new automaton rule where the prop- 
agation and velocity randomization steps are replaced by a mixing operator B that restores the particle 
distribution to a binomial distribution each automaton time step: 



Tr = C o 5 . 



(245) 



This automaton rule will allow us to explore the effects of fiuctuations on chaos in well-stirred systems 
where spatial pattern formation does not complicate the interpretation of the results. 



3. Simulation results 



The Willamowski-Rossler model has a chaotic attractor that arises by a period-doubling cascade. In 
all of the calculations presented below the rate constant parameters were fixed at the values: ki = 31.2, 
K_i = 0.2, K_2 = 0.1, K3 = 10.8, K_3 = 0.12, K4 = 1.02, K_4 = 0.01, K5 = 16.5, and k_5 = 0.5. The 
bifurcation diagram as a function of K2 for fixed values of the other parameters is shown in Fig. 19. The 
lattice-gas automaton has been used to explore the mesoscopic dynamics that underlies the subharmonic 
cascade leading to chaos. [66] 



67 



FIG. 19. Bifurcation diagram as a function of K2 showing period-doubling and reverse period-doubling cas- 
cades to chaos. The ordinate shows the concentration of species y in a Poincare section plane {pi = 
const. ,y p2 > 0,V/93 > 0}. The results were obtained by solving the mass action rate law (235). 



Well-stirred system 

Internal noise can have important effects on the dynamics when the system lies in or near the chaotic 
regime. From earlier investigations on the influence of external noise on nonlinear dynamical systems 
[69,70] one knows that noise can alter the positions of bifurcation points, destroy bifurcation sequences 
because periodic orbits cannot be resolved, induce periodic behavior and destroy the structure of chaotic 
attractors. All of these effects depend on the region in parameter space where the system lies and the 
amplitude of the noise and its statistical properties. Internal noise is somewhat more subtle to consider 
since it is generated by the system itself: the same dynamics that gives rise to the periodic or chaotic 
behavior on the macroscopic level is also responsible for the fluctuations. Therefore, the noise amplitude 
is not fully under control of the investigator and the possibility exits that the same mechanisms that are 
responsible for deterministic chaos may lead to an enhancement of the effects of fluctuations. It has been 
suggested that such effects can lead to the breakdown of the mean field macroscopic descriptions in the 
chaotic regime. [71] This is a rather paradoxical situation since the very equations that are typically used 
to analyse deterministic chaos may not be valid when the system is chaotic. The extent of this breakdown 
is a matter of debate. [72,73] 

The reactive lattice gas automaton can be used to explore these questions in a rather clear fashion. 
The automaton dynamics, while idealized, is based on the reactive dynamics as embodied in the reaction 
mechanism. Also, the mean field limit of the automaton dynamics is the mass action rate law. There- 
fore, the full automaton simulations which include fiuctuations can be compared with the mean field 
approximation to test its validity and the extent of its breakdown. 

It is especially interesting to study the dynamics for parameter values in the period doubling regime. 
Here one can see the effects of noise on limit cycle attractors of increasing complexity as one moves closer 
to and into the chaotic regime. Studies of this type have been carried out in Refs. [66] and [74]. We noted 
above that the fiuctuations scale as M'^^"^ but their actual magnitude depends on the system dynamics 
for the given control parameters. A series of automaton calculations has been carried out for parameters 
lying in the period doubling cascade for different system sizes. The two examples shown in Figs. 20 and 
21 illustrate the qualitative features of the internal noise effects. The period-two case is shown in Fig. 20. 
In this figure the period-2 attractors obtained from the automaton simulation for two different system 
sizes (middle and right panels) are compared with the deterministic period-2 orbit (left panel). 



FIG. 20. Projection of the 3-d phase-space trajectories on the (pi ,/92)-plane for K2 = 1.48 for the deterministic 
system (left panel) and for well-stirred automaton dynamics with A/" = (512)^ (middle panel) and A/" = (128)^ 
(right panel). Concentrations were obtained by averaging over all the nodes on the lattice. 

For small system sizes, below about M = (300)^, the noisy attractor bears little resemblance to the 
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deterministic period-2 orbit (cf. right panel). For large system sizes (middle panel) the noisy attractor 
is simply a thick version of the deterministic period-2 attractor in the left panel. 

The effects of fluctuations on the chaotic attractors are considered in Fig. 21 where the deterministic 
strange attractor is compared with the well-stirred automaton dynamics for two different systems sizes. 
For small system sizes = (1024)^, top panel) one observes that some of the fine structure of the 
chaotic attractor has been destroyed by the noise: all but the last the chaotic bands have merged and the 
dynamics has lost some of its phase coherence. In the middle panel results for M = (4096)^ are shown 
where additional fine structure of the attractor is resolved. 



FIG. 21. Chaotic attractors for two different system sizes: M = (1024)^ (upper panel) and A/" = (4096)^ (middle 
panel). Deterministic system (lower panel). P indicates Poincare section plane. 

In spite of these marked effects on the strange attractor, the gross structure of the chaotic phase space 
fiow is preserved. This is shown in Fig. 22 where the Poincare maps are displayed for the deterministic 
attractor along with those for two different system sizes. The noisy chaotic attractor has a line-like 
structure but is not banded for very small system sizes (A/" = (128)^, right panel of Fig. 22). 
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FIG. 22. Poincare maps showing the well-stirred automaton dynamics in the chaotic region, K2 = 1.568, for 
J\f = (128)^ (right panel), J\f = (1024)^ (middle panel) and deterministic system (left panel). 

Since the bands are no longer distinct, these map regions are no longer visited in a definite order on 
successive intersections of the phase space trajectory with the Poincare plane. As the system size increases 
a critical size is reached where the internal noise becomes sufficiently small that the chaotic bands are 
resolved and the noisy chaotic attractor appears to be a slightly perturbed version of the deterministic 
chaotic attractor (middle panel of Fig. 22). Only the structure on the smallest scales of the attractor 
hierarchy are destroyed by the noise. 

The noise scaling properties in the period doubling regime have been studied extensively for one- 
dimensional maps where a renormalization group treatment is possible. [69,75] Given equivalent points 
in the period doubling bifurcation sequence, for example the superstable points for period 2", one can 
ask how much the noise amplitude needs to be reduced to observe one additional period doubling. From 
numerical and renormalization group studies one finds that the noise amplitude must be reduced by about 
a factor of 6.7. Extensive calculations using the well-stirred automaton dynamics have verified that the 
number of particles must be increased by the square of this factor. So a similar scaling applies to both 
external and internal noise. [74] These results suggest that for macroscopic systems containing a mole of 
particles one should be able to resolve up to period 256 or 512 if no other sources of noise were present. 
This implies that in most chemical experiments external noise destroys the period-doubling structure 
before internal noise effects become important. However, for small systems this may not be the case. 
In addition, effects arising from spatial inhomogeneities can also lead to situations where internal noise 
effects can play a greater role. This will be discussed in more detail below. 

Spatially distributed system 

In an unstirred system diffusion is the only mechanism for removing local inhomogeneities in the chemical 
concentrations. If the system size is small enough and diffusion is strong enough, the spatial distribution 
of the concentration field will be nearly uniform and the results will be identical to those for the well- 
stirred system. If these conditions are not met, diffusion will be unable to homogenize the concentration 
distribution and local rather than global concentration fiuctuations will be most relevant for the dynamics. 
It is just such local fiuctuations that are responsible for nucleation-induced pattern formation processes. 

The results of a full automaton simulation of a small = (100)^) system are shown in Fig. 23 
(left panel) for K2 = 1.572, corresponding to the deterministic chaotic attractor in the right panel. 
(All simulations were carried out on triangular lattices.) Diffusion is able to smooth the concentration 
fiuctuations arising from reactions occurring on the time scale Tchem over a length L = (DTchemY^^ ■ For 
this automaton simulation L ~ 43, which is comparable to N = 100; thus, the system will maintain 
a roughly homogeneous spatial concentration distribution. This homogeneity was confirmed by direct 
observation of the distribution of chemical species over the lattices as a function of time. The full 
automaton dynamics is equivalent to well-stirred dynamics for this case. One observes in Fig. 23 that 
the noisy attractor has a similar gross structure to the corresponding deterministic chaotic attractor but 
differs in a number of respects. Internal fiuctuations cause the phase space trajectories to explore a 
greater volume of phase space leading to a larger size for the attractor. Furthermore, since the system 
size is small, fiuctuations are sufficiently strong to cause merging of the chaotic bands as noted above for 
well-stirred dynamics in small systems. Nevertheless, the density is non-uniform on the attractor so that 
regions of high probability density correspond to the underlying deterministic bands. Again this band 
merging is accompanied by dephasing of the dynamics when viewed in the Poincare plane so that the 
bands (or high density regions) are no longer visited in the same order. 
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FIG. 23. Chaotic attractor. Parameters are: ki = 31.2, k_i = 0.2, K2 = 1.572, k_2 = 0.1, K3 = 10.8, 
K_3 = 0.12, K4 = 1.02, K_4 = 0.01, K5 = 16.5, K_5 = 0.5. The dashed box which contains the attractor has one 
corner at the origin of the coordinate system, (x, y, z) = ( — .02599, —.05498, —.08193), and corners lying along the 
X, y and z axes at 2.51039, 2.25439 and 4.77463, respectively. The Poincare section is labeled P. Automaton 
simulation (left panel) and deterministic system (right panel). 

If the system is larger or the diffusion coefficient is smaller, diffusion may be insufficient to maintain 
spatial homogeneity over all of space. The simplest manifestation of such lack of diffusive mixing is the 
desynchronization of the chaotic oscillations in different spatial regions: local volumes (areas) of space 
containing many nodes oscillate in phase but the phase differs from local volume to volume. This is 
the analog of phase turbulence [76] for a periodic state but now the underlying oscillation is aperiodic. 
This effect is shown in Fig. 24. The left panel shows the chaotic phase space trajectory projected in 
the (pi , /92)-plane. As usual the concentrations Prik) were obtained by averaging over all nodes on the 
species lattices at each time k. However, the spatial distribution of particles is highly non-uniform as can 
be seen in the right panel. As a result, the trajectory occupies a small volume in phase space since the 
average over the lattice (real space) implies an average over different phases of the aperiodic attractor. 
In cases where spatial homogeneity no longer obtains, the effects of internal fluctuations are even more 
subtle since both diffusion and reaction combine to influence the magnitude of the local fluctuations. 



FIG. 24. Projection of the phase space trajectory in the (/9i/92)-pIane for the automaton dynamics with 
A/" = (100)^ and D = 1/40 (left panel). Spatial distribution of the Z = X3 species concentration is coded 
as gray levels (right panel). 

Provided the size of the well-stirred system is large enough, the noisy attractor, even in the chaotic 
regime, closely resembles that obtained from the deterministic rate law: only structure on the smallest 
scales is obscured by the noise. If the chaotic attractor has a hierarchical banded structure, as is the 
case for the Willamowski-Rossler attractor, there are critical values of the system size beyond which the 
structure at a given level in the hierarchy cannot be resolved. If the system size is small fluctuations can 
be so large that not only is the chaotic attractor signiflcantly modifled (cf. Fig. 23 - note that even in 
this case the gross geometrical structure of the attractor is preserved) but also the bifurcation sequence 
and locations of the bifurcation points can change. For example, noise can lead to premature truncation 
of the period-doubling cascade. 

If mixing is imperfect then spatial fluctuations in the concentration held also influence the dynamics. 
One can imagine the system to be composed of roughly homogeneous patches interacting with each other 
in a time-dependent fashion. Since the patches have linear dimensions much smaller than those of the 
entire system, the fluctuations within the patches have larger amplitudes. This feature, combined with 
the complicated but weak interactions among patches, leads to fluctuation effects that are strong and 
difficult to understand from simple considerations. 
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In experiments on macroscopic systems, the number of particles is so large that the effects of global 
concentration fluctuations are likely to be negligibly small, even in or near the chaotic regime where 
fluctuations are enhanced by the unstable dynamics. Since any experiment naturally involves some level 
of external noise, it would be difficult to distinguish internal and external noise effects: both act in a 
similar manner on the dynamics. [66] However, all experiments involve spatial degrees of freedom and 
local inhomogeneities. If diffusion is weak these spatial degrees of freedom come into play and the effects 
of fluctuations can become pronounced. In this regime one has to consider both the diffusion length scale 
as well as the temporal character of the dynamics. Thus, noise can influence the local structures that 
form in spatially-distributed systems. 
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X. MODELS WITHOUT EXCLUSION 



The exclusion principle does impose some limitations on the types of reacting systems that can be 
treated. For instance, there are limitations on the values of the rate constants that may preclude explo- 
ration of certain dynamical regimes. The origin of these restrictions was described in Sec. V where the 
scheme for the construction of the reaction probability matrix was given. One way to extend the range 
of applicability of the class of reactive lattice gas models is to modify or abandon the exclusion principle. 
Naturally there are penalties if this is done and some of the computational elegance is lost but systems 
can be considered that would be difficult to treat otherwise. There are a number of ways to construct 
reactive LGA models of this type but we focus on models that retain an integer- valued representation of 
the species concentrations [77-79] instead of real-valued representations [80] which lead to instabilities 
and are more akin to other finite-difference methods or coupled map lattices. Models using a Boolean 
representation but with a modified exclusion principle will not be considered here [31]. 

The reactive LGA models without exclusion studied thus far are similar. Basically particles reside 
at the nodes of a lattice with no restriction on the numbers of particles per node. Particles diffuse by 
hopping to neighboring nodes on the lattice and react according to probabilities that depend on the 
number of particles of a given species at the node. Since the reaction probabilities depend on the particle 
numbers and these numbers may in principle increase indefinitely, it is possible that the requirement that 
the reaction probability be less than or equal to unity will be violated. Conditions of the simulation can 
always be arranged so that this is an extremely rare event and is not encountered in finite space and 
time simulations. This restriction is a weaker version of that for the model with exclusion where it is 
appropriate to scale the dynamics so that the species concentrations rarely reach values corresponding 
to maximum occupancy of a node where exclusion distorts the transition probabilities and hence the 
fiuctuations. 

Briefiy the automaton rules may be summarized as follows. The rule used by Chopard et al. [77] may be 
composed of diffusion and reaction steps. Diffusion is simulated by allowing particles to jump to nearest 
neighbor sites with equal probabilities for each lattice direction. The diffusion coefficient is varied by the 
number Ij of consecutive diffusion steps for a given species r so that Dj = Ij /'^ (for a square lattice). 
Reactions of a specific type are considered. For a reaction of the form niXi -\- U'^X'^ 113X3, all possible 
reactions composed from ni Xi and n2 X2 particles have a probability k to create X3 particles. The 
automaton rule considered by Karapiperis and Blankleider [78] is similar and again consists of diffusion 
and reaction steps. Diffusion is simulated by allowing particles of species r to move to neighboring lattice 
nodes with a given probability. The probabilities may be different for the various species and are chosen 
to satisfy isotropy conditions or may be biased depending on the application. Reactions are carried out 
with probabilities that depend on the numbers of particles of the given species at a node. The "rule IF' 
in [78] is similar in form to that described earlier for reactive LGA and leads to the standard mass action 
rate law when correlations are neglected. 

We illustrate the application of such reactive LGA models without exclusion using the following rule 
[79]: The dynamics is carried out on a stack of species lattices with coordination number m^. The number 
of channels at each node is taken to be equal to the coordination number (m = m^; allowed velocities 
correspond to jumps from node to nearest neighbors node in a discrete time step). Now there is no 
restriction on the numbers of particles that can occupy a given channel. 

• Diffusion is simulated by deterministic motion of particles from a node to its neighbors along 
the directions specified by their velocities. Following propagation the velocities are randomized 
by reassigining the particles with equal probabilities to the m velocity directions. As noted above, 
there is no restriction on the number of particles with a given velocity. These are just the analogs of 
the propgation and velocity randomization processes described earlier for the model with exclusion. 

• Reaction is carried out by probabilistic particle number changes similar to those described for the 
model with exclusion except that the probability distributions are Poisson instead of binomial. The 
elements of the probability matrix can be determined from the formulas given earlier by replacing 
ratios of factors involving m by unity, i.e., in the elements of P{oL,(j) the limit m ^ cxd is taken. 
As in the abovementioned two models, it is possible in the course of simulation that the computed 
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reaction probability will exceed unity but this can always be arranged to occur with negligible 
probability and when it does occur the corresponding probability is set equal to zero. 

As an illustration of the application of the reactive LGA without exclusion we consider the Willamowski- 
Rossler model described earlier. [67] Figure 25 shows two noisy limit cycle simulations carried out on small 
lattices (100 x 100). Panel (a) is the result using an LGA model with exclusion while panel (b) shows 
the result for the automaton without exclusion. 



FIG. 25. Comparison of noisy limit cycles for the WR model, (a) with exclusion and (b) without exclusion. 
The parameters are: K2 = 1.4, D = 1, h = 1.5 X 10-\ The origins of the boxes containing the attractors are at 
(0.10440,0.02660,0.2600) while the opposite corners lie at (1.24540,1.55900,3.86590) 

The fact that the two simulation methods give similar results is not surprising in view of the fact that 
the average particle number densities are such that maximum occupancy of a node is a rare event. In 
fact, if in the models with exclusion this is not the case the fluctuations predicted by the model should 
be suspect. As noted earlier, the model without exclusion allows a larger class of chemical systems to 
be treated with a corresponding increase of computational complexity. Thus with minor modifications 
the methods and algorithms developed in the preceding sections can be carried over to models without 
exclusion should the need arise. 
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XI. CELLULAR AUTOMATA 



The reactive lattice-gas automaton is only one of a class of discrete models that may be constructed 
for the description of spatially-distributed chemically reacting systems. Such discrete models include, 
among others, cellular automata (CA), coupled map lattices, lattice Boltzmann equations and birth- 
death master equation models. We shall not consider coupled map lattices further since the dynamical 
variables are continuous; comments on their relation to CA and reactive LGA can be found in Ref. [50] 
and their application to reactive systems is described in Ref. [81]. Lattice Boltzmann methods [82] also 
utilize real- valued dynamical variables and in this sense lie in the same class of systems as coupled map 
lattices. Their application to reactive systems can be found in [83]. As noted in the Introduction and 
in Secs.VI-IX, master equation methods are closely related to the lattice gas methods described here; 
discussions of their application to reacting systems are given in Refs. [11-13] as well in the specialized 
references quoted in Secs.VI-IX. We shall confine our attention to cellular automaton models since these 
are closest in structure to the LGA method. In this section we contrast the lattice gas automaton scheme 
with more traditional implementations of cellular automaton ideas. 

Cellular automata are abstract discrete dynamical systems devised by Von Neumann to describe evo- 
lution in biological systems. [84] A cellular automaton consists of a set of nodes, usually arranged on a 
regular lattice; each node supports state variables that take on a finite number of possible values. The 
state variables are synchronously updated at discrete time intervals according to a local deterministic or 
probabilistic rule that depends on the state variables at and in a neighborhood of a node. It is clear that 
reactive LGA lie within the general class of models given by this prescription and can formally be called 
automata (for a different viewpoint see [85]); however, the modeling strategy and rule construction is 
quite distinct from that of traditional cellular automaton models since the aim of reactive LGA is the 
construction of a mesoscopic description of the reactive collision dynamics and the rule is composed of 
propagation and collision steps. 

It is now well established that the above simple prescription for the rule construction can lead to CA 
with a bewildering variety of complex dynamics, even for Boolean state variables and a one-dimensional 
lattice of nodes with nearest-neighbor interactions. There exists a large literature on the formal properties 
of such automata. [86] Cellular automata can be constructed as simplified models of reaction-diffusion 
systems. Often the main features of an apparently very complex dynamics can be captured in a very 
simple rule. Cellular automaton models have been constructed for a wide variety of chemically reacting 
systems. This section is not intended to be an exhaustive review of the extensive literature in this field; 
instead a few examples will be given to illustrate and contrast the point of view taken in the CA model 
construction with that of reactive LGA. Further references to the chemical CA literature can be found 
in Ref. [50]. 

Often the CA modeling strategy can be described as the inverse problem: given a phenomenon, what is 
the cellular automaton rule that produces it? [87] One generally starts from some notion of the physical 
ingredients that characterize the system and builds them into the rule. A good rule will yield dynamics 
that mimics that of the real system. Once a class of rules that depend on parameters, such as the 
neighborhood size or the number of states, is constructed, it is interesting to attempt to classify the 
possible types of dynamical behavior that arise from such parameter variations. This is the forward 
problem in cellular automata theory. Since cellular automaton rules are notorious for their lack of 
smooth behavior under such parameter variations, it is a topic of considerable interest in applications. 

Consider, for example, excitable media discussed in Secs.VI-IX. The reactive lattice gas modeling of 
such excitable media is really no different from that for any other type of reactive dynamics: one starts 
with the reaction mechanism and constructs a microdynamics that mimics the reactive and nonreactive 
collision processes in the system. If the kinetics happens to give rise to excitability then so will a properly 
constructed reactive lattice gas model. This was illustrated in Secs.VI-IX where the Selkov model was 
studied. The same general Selkov lattice gas microdynamics was shown to describe excitability, bistability, 
Turing pattern formation, etc. Thus, no special features of excitability are exploited in order to construct 
the automaton rules. The situation is different in traditional applications of cellular automata and below 
we outline some aspects of the way one constructs CA models for such media. 

Excitable media were among the first systems to be modeled using cellular automaton methods. In 
an early paper Wiener and Rosenblueth [88] constructed a cellular automaton-like model for excitable 
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cardiac tissue. They recognized the principle features of an excitable medium described in Secs.VI-IX: 
an excited state at the wave front, refractory states where the medium can "neither transmit or receive 
impulses" and a resting state that follows the termination of the refractory state where the system is 
susceptible to stimulation. While these basic ideas and their variants were used in subsequent studies of 
excitable media (cf. Farley's work on neural activity [89]), the most extensive early investigations of these 
models were carried out by Greenberg and Hastings [49]. The subsequent generalizations of their work 
have been the subject of a great deal of mathematical study and form the basis of recent applications. 
The Greenberg-Hastings (GH) rule construction illustrates how the above characterization of an excitable 
medium can be transcribed into a simple CA rule with rich dynamics, only some of which resembles that 
of real excitable media. 



A. Cellular automata for excitable media 

The Greenberg-Hastings model can be written as the following rule: the dynamical variable u describing 
the state of the system can take on the three values u = { — 1, 0, 1} where is the resting state, 1 is the 
excited state and —1 is the refractory state. The three-state excitable medium cellular automaton rule 
is: 

rO, ifM(r,A;) = -l 

u{r,k+l)= I -1, ifM(r,A;) = l . (246) 

I max{0, M(r', A;) : r'eW'(r)}, if M(r, A;) = 

Here r labels the sites of a regular lattice and Af(r) is a neighborhood of r. This very simple rule embodies 
the essential features of excitable medium dynamics. In the context of this rule it is simple to understand 
the formation of the two principal patterns of excitation, rings and spiral waves, and the initial states 
that give rise to their formation. It is also possible to appreciate some of the complexities of the dynamics 
starting from random initial conditions in terms of this rule. 

Generalizations of this simple model are termed GH rules. Again, let r label the nodes of a regular 
lattice. The discrete dynamical variable u takes n integer values in the set {0, . . .,n — 1}. A site can 
change its value either independently of the values of its neighbors or in a manner that depends on the 
values of the neighboring sites. The interaction with neighboring sites is characterized by two additional 
parameters: a threshold 6 and a range 6. Let Afs(r) = {r' : |r' — r| < (5} where | ... | is a suitably defined 
norm. Let Ns be the cardinality of the set {r' G Afs(r) : u(r' , k) = (u(r,k) + l)mod n}. With these 
definitions the GH rules are 

M(r fc-Fl) = 1'^"'^^'^^ + ^^™°^"' if ^) > 1 or #i > 6* .^47) 
^ ' ' ^u(r,k), otherwise ' ^ ' 

This is the generalization of the simple n = 3, 6 = I, 6 = 1 rule originally studied by GH. This generalized 
GH rule has been the subject of extensive mathematical investigations. [90-92] Some aspects of its ergodic 
properties are known and the broad structure of the phase diagram as a function of n, 6 and 6 has been 
determined. Fisch, Gravner and Griffeath have shown that GH rules can exhibit a variety of asymptotic 
states depending on the rule parameters and, in addition, the relationship between these models and 
continuum reaction-diffusion models for excitable media has been examined. [90] Probabilistic versions 
of excitable medium CA rules have also been constructed and investigated. [90-92] In addition to these 
simple GH rules, more complex rules that are designed to model specific features of the wave propagation 
processes in excitable media have been constructed. In particular, we mention the work of Gerhardt, 
Schuster and Tyson [52] and Markus and Hess. [51]. These CA models can account for the dispersion of 
the wave velocity and, in the case of the latter model, the wavefront curvature. 

B. Cellular automata for reaction-diffusion systems 

CA for reaction-diffusion systems are discrete models which offer an alternative to partial differential 
equations. They have been shown to produce qualitatively correct behavior. However, they are often 
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restricted to certain RD systems and certain types of phenomena and have been subjected to the criticism 
of experimentalists and researchers working with PDEs who are concerned with quantitative predictions. 
A class of CA was developed for modelling many reaction-diffusion systems in a quantitative way [93]. 
The main idea behind this class of CA is careful discretization. Space and time are discretized as in 
normal finite difference methods for solving the PDEs. Finite difference methods then proceed to solve 
the resulting coupled system of N x n ordinary differential equations (N points in space, n equations 
in the PDE system) by any of a number of numerical methods operating on fioating point numbers. 
The use of fioating point numbers on computers implies a discretization of the continuous variables. The 
errors introduced by this discretization and the ensuing roundoff errors are often not considered explicitly, 
but assumed to be small because the precision is rather high (8 decimal digits for usual fioating point 
numbers). In contrast, in the CA approach, all variables are explicitly discretized into relatively small 
integers. This discretization allows the use of lookup tables to replace the evaluation of the nonlinear rate 
functions. It is this table lookup, combined with the fact that all calculations are performed using integers 
instead of fioating point variables, that accounts for an improvement in speed of orders of magnitude on 
a conventional multi-purpose computer. The undesirable effects of discretization are overcome by using 
probabilistic rules for the updating of the CA. 

The state of the CA is given by a regular array of concentration vectors p residing on a rf-dimensional 
lattice. Each p(r) is a n- vector of integers, where n is the number of reactive species and each component 
/9^(r) can only take integer values between and where the &^'s can be different for each species r. 
The position index r is a rf-dimensional vector in the CA lattice (for cubic lattices, r is a rf-vector of 
integers). 

The central operation of the automaton consists of calculating the sum of the concentrations in some 
neighborhood Afr 

Pr{r)= Pr(r + r'). (248) 

The neighborhood (which can be specific to each species r) is specified as a set of displacement vectors, 
e.g. for a two-dimensional square lattice the neighborhood restricted to first neighbors is given by 



For some - in particular square and cubic - neighborhoods, the summation operation can be executed 
very efficiently [93]. It is convenient to introduce the normalized values £'r(i') = Pr(i')/&r and ^r(i') = 
p^(r)/(&^|jVV|), which are always between zero and one. The resulting fields ^r(i') are then the local 
averages of the £>r(r). The averaging has the effect of diffusion as can be seen from a Taylor expansion 
of gr(r + r') around £>r(r): 



' ^ ' r'GA^ ;=0 ^ 

= er(r) + i?.V2e.(r) + (250) 

The factors Dj can be computed as in [94] and are easily calculated from (250) for square neighborhoods 
with radius I: Dr = 1(1 + l)/6. 

The second operation in the cellular automaton is the implementation of the reactive processes described 
by a rate law. Given the reaction-diffusion equation 

^ = F{p) + VV'p, (251) 

we discretize the time derivative to obtain 

p*+^* = p* + AtF{p*) + AtVV-'p*. (252) 
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Changing the time and space scales by setting k = t/ At and r = x/ Ax, and using the variable g for the 
rescaled set yields 

e''+' = + AtFiQ'') + -p-VV'QK (253) 

(Ax)^ 

as the equation to be treated by the CA. 
Let us define 

F*{q) = Q + AtF{Q). (254) 

From (250) and (254) 

F*{q'') = q^ + T>V'q^ + • • • + AtFie^ + T>V'q^ + • • •) 

= +Y)V'^Q^ + AtF{Q^)+ OiAt"^). (255) 

Then 

q^+^ = F*{q^) (256) 

is consistently first order accurate in time and, within this limit, (256) can be validly identified with (253) 
to describe the evolution of the system. The identification yields 

Dr = ^Vr (r= !,•••,«) (257) 
Ax 

which defines the space scale. As p is the result of the diffusion step, the average output of the CA 
reaction-diffusion process is therefore given by 

for species k. Probabilistic rules are important. Given an input configuration p*^(r), one assigns new 
values p^'^^{v) probabilistically in such a way that the average result corresponds to the finite difference 
approximation to the given reaction-diffusion equation, Vj (for further details see [95]). 

An interesting application is the mapping of the Ginzburg-Landau equation onto the automaton. Con- 
sider the PDE [76] 

^ = DV"^ z + az - h\z\'^ z , (259) 

which can be viewed (for D real) as a two-species reaction-diffusion system by separating real and 
imaginary parts z = x -\- iy , a = a -\- ij, and h = [3 -\- ib, 

^ = DV'^x + ax-jy + {-px + 8y){x'^ + t/^) , (260) 
^ = DV^y + ay + jx + {-py - bx^x-" + j/^) , (261) 

where x, y are space- and time-dependent variables. With z = re"^, the corresponding rate equation can 
be conveniently written as 

= ar - (]r^ = (3r {a/ (3 - r^) , (262) 
^-^ = l-8r\ (263) 
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One can set /3 = 1 by changing the time scale. The stable homogeneous solution is z = rje'^* with 
= \/ aj fi and SI = 7 — bri (A steady but unstable solution is z = 0). Since the equation for r is 
independent of one can transform the solution for SI = to any given by multiplying z{t) with 
a factor e"^*, yielding an oscillating solution. The full PDE system, (259), also admits oscillating and 
rotating spirals, and other inhomogeneous solutions. 

Study of the homogeneous solutions shows that the CA behaves like an explicit finite difference method 
with added noise. In the automaton, the noise is intrinsic (it arises from the discretization) but has little 
effect on the amplitude of the solution (although it introduces random drifts in the phase). An interesting 
situation concerns spiral wave solutions. By starting with the initial condition 

i?(r,, r,, t°) = ci^(r, - rlf + (r, - r|;)2 , (264) 

6'(ra;, ry,t°) = C2 arctan ^_ g , (265) 
I'x — 

which creates exactly one phase singularity at (r^, rJJ), for nonzero 6 smaller than a critical value, a spiral 
develops and rotates steadily after some time as shown in Fig. 26(left panel). The value of the wavelength 
(obtained by an Archimedian spiral fit) is in agreement with the theoretical value [96] ( for small enough 
At). 



FIG. 26. Spiral waves in the Ginzburg-Landau system obtained by CA simulation. Left panel: Automaton size 
200^ sites (Superimposed on the grayscale plot of x is the fitted Archimedian spiral). Right panel: Automaton 
size 500^ sites. 

As an example of greater complexity. Fig. 26(right panel) shows the simulation of a large system (500^ 
sites) initialized in the state z 0: Under such conditions many interacting spirals develop. Indeed, the 
initial state is unstable, and because of the intrinsic (low level) noise, different regions depart from this 
unstable state with different phase values (f). This situation automatically creates many phase singularities 
(points with r = 0, surrounded by points with all values of (f)), which then develop into spirals. These 
phase singularities can merge and move as they are infiuenced by each other [97]. Such large system 
simulations are made possible by the speed advantage that this class of CA - as compared to the class of 
LGA 2^ - offers over other numerical methods for solving PDEs. Further applications include the study 
of spiral patterns [98] and the spatial coexistence of different patterns [99]. 



^ Note however that LGA simulations performed on the CAM-8 machine compare in computation speed to PDE 
solvers. 
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C. Remarks 



The brief discussion in this section was intended to illustrate how the reactive lattice gas method differs 
from both traditional and more modern CA modeling of reaction-diffusion systems. In the reactive LGA 
method the dynamics is modeled at the mesoscopic level where one attempts to construct a fictitious 
albeit faithful dynamics that describes the reactive and non-reactive collision events. The reaction prob- 
ability matrix that controls the reactive dynamics may be built on the reaction mechanism and suitably 
constructed to yield the reaction-diffusion equation in the mean field limit. If this matrix is properly 
constructed as described in Sec. V, then particle number fiuctuations that arise from reaction and dif- 
fusion will be described correctly. Features such as wave dispersion and curvature are automatically 
incorporated and follow from the mesoscopic dynamics. No special additions to the rule are needed to 
achieve these results. 

Of course, since the model building takes place at the level of the mechanism the reactive LGA rule is 
as complex as the mechanism. If one is interested solely in the macroscopic wave propagation properties 
then a PDE level of description will suffice and CA models, even rather complex ones, may provide 
a computationally efficient way to simulate this macroscopic reactive dynamics. The results in the 
preceding subsection describe one such method. However, as stressed elsewhere in this paper, the reactive 
LGA allows one to treat fiuctuations in a realistic manner and can be used to explore regimes that are 
inaccessible by PDE or simple CA methods. 
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XII. PERSPECTIVES 



There are circumstances when only the global behavior of spatially-distributed reacting systems on long 
time and distance scales is of interest. In this domain macroscopic reaction-diffusion equations provide an 
appropriate way to model the dynamics and, provided the systems are sufficiently simple, the boundary 
value problems implied in the solutions of these partial differential equations can be implemented and 
solved numerically. However such a macroscopic model is not capable of describing the system on short 
distance and time scales, nor can it account for microscopic fluctuations which can be amplifled by 
the dynamics to macroscopic scale. Since reaction-difusion equations are "mean fleld" in character they 
cannot give a full description of the correlations which can lead to the breakdown of standard mass action 
rate laws. Also, the systems of interest may have complex, inhomogeneous structure, as in porous media, 
and it may be impossible or very difficult to implement the boundary conditions on the reaction-diffusion 
equations in such cases. 

The reactive lattice gas automata described in this paper are designed to provide a microscopic approach 
to complex, spatially-distributed, reacting systems and allow an analysis of these systems at a mesoscopic 
level. Indeed it is an ideal method for the exploration of reactive dynamics on the interesting scales that 
lie between the microscopic, where full molecular dynamics must be used, and the macroscopic, where 
reaction-diffusion equations suffice. In addition, it can be used to explore the dynamics of systems near 
bifurcation points or in chaotic regimes where fluctuations can have important effects. 

Since the method is based on a well deflned reactive particle dynamics with integer-valued state vari- 
ables, it does not suffer from any of the instabilities of traditional flnite-difference schemes. It is easily 
applied to systems with complex boundaries since one need only include collisions with the "walls" to 
account for their presence. [101] Since the method is akin to a molecular dynamics simulation in that 
particle dynamics is followed, simulations of a given system (with simple geometry) are generally slower 
than PDE or traditional cellular automaton models - with the advantage that LGA provide a much 
deeper level of description of the system. Nevertheless the reactive LGA can be implemented very effi- 
ciently on parallel machines as well as special purpose machines such as the CAM-8 [102] and in these 
circumstances the speed of LGA simulations may surpass that of reaction-diffusion PDEs. However, we 
stress that the real utility of these reactive LGA methods lies in the more fundamental treatment of the 
system that they provide. 

The applications in sections 6 to 9 show the types of new information that can be obtained from the 
study of reactive LGA models. In particular we should stress a few important points: 

- Simulation results and analytical developments conflrm the validity of the phenomenological descrip- 
tion of Reaction-Diffusion systems when the Boltzmann hypothesis is applicable; when the frequency of 
reactive collisions becomes large, non-Boltzmann effects are observed and the phenomenological descrip- 
tion breaks down; 

- LGA provide qualitative insight in nucleation processes and in the early stages of pattern formation; 

- the role of fluctuations can be studied quantitatively; 

- LGA results on reactive chaos provide a physically motivated basis to phenomenological treatments; 

- reactive LGA have the same level of validity and fundamental character as Master Equation models 
with the advantage that with LGA we perform mesoscopic simulations on large, spatially-distributed sys- 
tems which could hardly be realized with Master Equation methods and Molecular Dynamics techniques. 

Many problems which could be investigated from the LGA approach remain unexplored. For instance, 
very little work has been done on the application of these methods to reactive flows ^® and extensions 
of the theory are necessary before this problem can be fully mastered. Another question which remains 
open is the incorporation of energy levels in the automaton - as in thermal LGA [40] - to account for 
the kinetic processes at a foundamental level. The class of systems that can be treated by the methods 
described here is in fact much broader than the speciflc applications indicated. An interesting extension is 
the modeling of polymerization through hetrogeneous catalysis. [101] Many diverse problems can actually 



Recent work on nonlinear reactions advected by a flow using the Lattice Boltzmann method shows how flows 
can modify the resulting effects of the reactive processes. [100] 
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be cast into the form of a reactive dynamics. For example, the dynamics of the populations of excited 
atomic or molecular states in lasers can, in some circumstances, be modeled as reactions between different 
"chemical" species. [103] There is an even broader class of phenomena in biology that lend themselves 
to such modeling. One may generalize the notion of automaton particles to have them represent, for 
example, motile cells (bacteria or amoebae [104,105]) undergoing "random walks" yet showing collective 
behavior at the macroscopic level, just as the molecules in some of the chemical examples treated in this 
paper. Thus, the reactive lattice gas automaton can provide a theoretical and computational method for 
the study of diverse classes of phenomena at the mesoscopic scale. In some of these classes there are open 
questions where the automaton approach offers interesting perspectives. 
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